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ABSTRACT 


This  thesis  investigates  the  applicability  of  VRTs 
to  the  simulation  of  stochastic  combat  models.  Ways  of 
measuring  the  efficiency  of  a  VRT  are  explored. 
Antithetic  variates  and  stratified  sampling  are  applied 
to  the  simulation  of  a  trivariate  Markovian  combat 
model.  Means  of  programming  the  antithetic  variates  a.  c 
stratified  sampling  to  reduce  the  inherent  variability 
of  uncertainty  in  the  output  data  of  the  model  are 
illustrated.  Response  surface  regression  models  are 
used  to  characterize  the  performance  of  the  antithetic 
variates  and  stratified  sampling  in  the  Markovian 
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combat  model. 
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INTRODUCTION 


I  . 


A .  GENERAL 

War  policies  and  plans  for  military  operations  made 
during  peacetime  are  significant  for  the  mission 
accomplishments  of  combat  operations  conducted  during 
wartime.  Since  experimentation  with  real  combat  is 
infeasible,  military  analysts  use  stochastic  combat 
simulation  models  to  study  the  effects  of  policy  making 
on  combat  operations.  The  analysts’  inferences  drawn 
from  the  results  of  these  models  are  important  to  the 
decision  maker  since  he  has  to  use  them  to  make  the 
same  decisions  about  military  operations  as  he  would  if 
he  could  experiment  with  real  combat  itself.  The  output 
data  from  these  models  are  realizations  of  random 
variables  distributed  around  the  values  of  the 
parameters  of  interest,  or  the  models’  true 
characteristics,  so  the  analysts  can  only  estimate 
these  parameters  with  error.  The  magnitude  of  error 
for  each  estimate  can  be  measured  in  terms  of  precision 
or  the  variance  of  the  estimate:  if  the  estimate  is 
unbiased  then  the  smaller  the  variance,  the  greater  the 
precision,  and  the  smaller  the  error.  Since  the 
decision  maker  wishes  to  make  decisions  that  are  based 
on  estimate(s)  with  a  quantifiable  error  bound,  the 
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analysts  may  find  it  possible,  to  apply  specific 
statistical  techniques  to  measure  and  control  the 
variance  of  the  estimate(s)  to  obtain  a  prescribed  or 
at  least  quantifiable  level  of  precision. 

The  analyst’s  capability  of  estimating  a  parameter 
of  interest  with  high  precision  depends  on  the  extent 
to  which  he  is  able  to  control  the  sample  variance. 
When  the  analyst  uses  the  mean  value  of  the  sample 
output  data  as  the  estimate  of  a  parameter  of  interest 
and  when  individual  samples  are  independent,  the 
coefficient  of  the  variance  of  the  estimator  is  reduced 
by  a  factor  of  1 /n  where  n  is  the  sample  size  of  the 
output  data.  A  large  sample  size  yields  an  estimate 
with  a  small  variance  and  high  precision.  Multiple 
replications  to  obtain  a  large  sample  size  in  complex 
stochastic  models  can  be  prohibitively  expensive  in 
terms  of  resources  like  money,  internal  computer  time, 
computer  storage  space,  etc.  This  is  especially  true 
for  large-scale,  complex  stochastic  combat  simulation 
models,  which  often  require  hours  rather  than  minutes 
for  a  single  computed  replication.  Since  available 
computer  time  is  a  compelling  constraint  on  military 
studies  competing  for  scarce  resources,  the  analyst  is 
usually  given  an  allocated  amount  of  time  to  simulate 
his  model.  This  specified  amount  of  time  may  affect  a 
desired  level  of  precision  of  the  estimate(s)  that  the 
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analyst  wishes  to  obtain  from  the  simulation.  Since 
the  analyst  can  only  execute  a  fixed  number  of 
replications  within  this  block  of  time,  the  sample  size 
(number  of  replications)  may  not  be  large  enough  to 
achieve  a  variance  small  enough  to  give  the  analyst  an 
acceptable  precision  for  the  estimate(s).  Hence,  the 
analyst  must  either  accept  the  particular  level  of 
precision  and  error  associated  with  such  variance  or 
apply  other  specific  statistical  techniques  which  are 
more  likely  to  produce  a  smaller  variance,  and  hence  a 
level  of  precision  with  which  he  can  feel  more 
comfortable . 

An  economizing  scheme  in  simulation  to  reduce  the 
variance  of  the  estimator  is  to  intentionally  distort, 
control,  and  modify  the  random  properties  of  the  input 
variables  in  the  simulation  model.  The  output  data 
resulting  from  the  manipulation  of  these  random  numbers 
are  random  variables  which  are  designed  to  be  much 
closer  together  and  more  closely  distributed  around  the 
true  value  of  the  model’s  parameter  of  interest  than  is 
the  case  with  simple  random  sampling.  A  sample 
distribution  resulting  from  such  a  variance-reducing 
scheme  has  the  same  mean  value  but  a  potentially 
smaller  variance  than  the  distribution  of  the  sample 
without  the  usage  of  this  scheme.  The  different 
techniques  for  doing  this  scheme  are  called  Variance 
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Reduction  Techniques  (VRTs).  The  effects  of  certain  of 
these,  when  applied  to  a  combat  model,  are  the  subject 
of  this  thesis. 

B.  BACKGROUND 

VRTs  were  initially  used  to  evaluate  multi¬ 
dimensional  integrals.  They  nave  since  been  applied  to 
small  Monte  Carlo  simulation  problems  but  have  not  been 
extensively  utilized  in  large  complex  stochastic 
simulation  models.  The  utilization  of  these  variance- 
reducing  techniques  in  real-wo’-ld  combat  simulation 
models  is  even  less  common  Consequently,  limited 
examples  of  applications  of  VRTs  in  these  simulation 
models  are  found  in  the  literature.  The  major  reason 
for  this  is  because  the  performance  of  the  VRTs  is 
suspected  to  be  uncertain  and  unpredictable.  The 
analyst  has  no  guarantee  that  the  usage  of  VRTs  will 
work  all  the  time.  Futhermore,  he  has  no  way  to  know 
beforehand  how  much  variance  reduction  he  will  get  from 
the  application  of  VRTs  whenever  they  are  effective. 
However,  VRTs,  in  our  opinion,  promise  to  be  powerful 
and  effective  tools  in  simulation  if  the  issues  of 
their  performance  in  specific  simulations  are 
understood.  In  this  section  we  will  describe  the  effect 
that  they  can  have  on  simulation  studies.  In  Chapter  V 
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we  illustrate  their  effectiveness  for  a  particular 
combat  model. 

The  effectiveness  of  a  VRT  may  be  measured  by  the 
relative  efficiency  of  the  simulation  in  obtaining 
estimate(s)  with  the  utilization  of  this  scheme,  to 
the  efficiency  of  a  simulation  under  the  same 
conditions  without  the  VRT.  Efficiency  as  Handscomb 
(1969,p.  253)  defines  it  is 

Efficiency  =  1  /  (Variance  *  Work).  (i ) 

Here  "Work"  generally  refers  to  computing  time. 
According  to  Handscomb,  variance  reduction  succeeds  if 
the  VRT  increases  efficiency.  From  Equation  1,  we  see 
that  a  decrease  in  variance  and/or  work  will  increase 
efficiency.  Hence,  variance  reduction  in  simulation  is 
more  than  solely  a  decrease  in  the  variance  of  the  mean 
of  the  estimators.  Handscomb( i 969 ,p .  253)  calls  a 

technique  var lance -reducing  if  it  "reduces  the 

variance  proportionately  more  than  it  increases  the 
work  involved"  or  "does  not  reduce  the  variance  at  all 
in  the  usual  sense,  provided  that  it  saves  work."  The 
work  involved  in  attaining  estimates  by  simulation  has 
many  attributes.  Hamraersley  and  Handscomb( 1 964 , p .  22) 
suggest  that  the  number  of  simulation  runs  epitomizes 
this  work.  However,  we  can  easily  measure  this  same 
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work  in  terms  of  computing  cost  or /and.  simulation  time. 
For  it  is  the  availability  of  these  factors  that 
ultimately  determine  the  precision  of  the  estimators. 
Hence,  an  effective  VRT  may  not  only  produce  more 
precise  estimates  but  also  economize  the  time  and  costs 
associated  with  the  simulation  to  obtain  the  level  of 
precision  for  those  estimates.  The  efficiency  of  VRTs 
will  be  discussed  more  fully  in  Chapter  III. 

This  potential  saving  in  computer  time  has 
stimulated  the  interest  of  the  United  States  Army 
Concept  Analysis  Agency  ( CAA )  in  the  utilization  of 
VRTs.  CAA  has  studied  the  effectiveness  of  a  VRT  in  two 
of  its  larger  and  more  complex  simulation  models.  The 
results  of  these  studies  were  mixed  (Johnson,  Bates, 
and  Graham,  1985).  CAA  therefore  recommended  the 
continuation  of  the  studies  to  investigate  the 
applicability  of  VRTs  to  reduce  the  inherent 
variability  in  large,  complex,  stochastic  combat 
simulation  models. 

C.  PURPOSE  AND  OBJECTIVE 

The  purpose  of  this  thesis  is  to  provide  additional 
insight  into  the  applicability  of  VRTs  to  stochastic 
combat  models,  and  to  provide  a  base  for  future  studies 
in  the  application  of  VRTs  to  large-scale,  real  world, 
stochastic  combat  simulations.  The  objectives  of  this 


thesis  are  plainly  to  identify  those  VRTs  that  are 
applicable,  and  then  to  exhibit  their  performance  in 
the  applications  to  a  class  of  stochastic  combat 
simulation  models.  The  question  to  be  answered  is:  "Can 
VRTs  be  identified  that  are  consistently  effective  for 
reducing  simulation  time  and  cost?" 

D.  PROBLEM 

The  problem  for  this  thesis  is  to  increase  the 
efficiency  of  a  stochastic  combat  simulation  model 
utilizing  VRTs  in  terms  of  (1 )  increased  precision  of 
the  model’s  estimates  for  an  allocated  amount  of 
simulation  time,  and  (2)  reduced  computer  time  for  a 
predetermined  level  of  precision. 

E.  APPROACH 

In  his  doctoral  dissertation,  Andrtasson  (1972) 
shewed  that  variance  reduction  in  queuing  systems  is 
influenced  by  (i)  the  transformation  of  random  numbers, 
(ii)  the  structure  and  parameters  of  the  simulation 
model,  and  (iii)  the  choice  of  the  model  response 
quantity.  Condition  (i)  is  an  attribute  of  the  VRTs. 
Conditions  (ii)  and  (iii)  are  characteristics  of  the 
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this  thesis,  on  variance  reduction.  We  then  use  those 
results  to  formulate  our  approach  to  increase  the 
efficiency  of  this  model  in  terms  stated  in  the  blem 
above . 

F.  ORGANIZATION  OF  THIS  THESIS 

This  thesis  is  organized  into  6  chapters.  Chapter  I 
is  the  Introduction  chapter.  Chapter  II  reviews  the 
literature  of  VRTs  in  simulation.  Chapter  III  discusses 
ways  of  measuring  the  efficiency  of  a  VRT  and  explores 
the  tradeoffs  of  measuring  for  increased  precision  of 
estimation  and  reduced  computer  time.  Chapter  IV  is 
concerned  with  the  simulation  of  a  stochastic  combat 
model  and  the  programming  for  variance  reduction  in  the 
simulation  model.  Chapter  V  deals  with  the 
applicability  and  performance  of  VRTs  in  the  simulation 
model.  In  Chapter  VI,  we  make  conclusions  about  the 
applicability  of  VRTs  in  stochastic  combat  models  and 
provide  recommendations  about  their  use  in  larger  and 
more  complex,  stochastic  models  that  are  used  to  study 
real-world  combat  systems. 


II .  REVIEW  OF  LITERATURE 


A.  INTRODUCTION 

The  VRTs  that  we  use  to  solve  the  problem  stated 
In  Chapter  I  of  this  thesis  are  antithetic  variates  and 
stratified  sampling.  But  first  we  review  the  literature 
of  variance  reduction  ir.  simulation.  This  chapter 
concentrates  on  the  practical  applications  of  VRTs  in 
simulation  models.  We  present  a  brief  summary  of  works 
of  scholars  and  experts  on  this  subject.  We  then 
describe  the  basic  concepts  of  two  VRTs  that  we  feel 
are  applicable  to  large-scale,  complex,  stochastic 
combat  simulation  models.  It  is  these  two  VRTs  whose 
performances  we  later  exhibit  in  the  combat  model  in 
this  thesis. 

B.  SUMMARY  OF  PREVIOUS  WORKS 

In  the  last  15  years  Interest  in  VRTs  In  simulation 
has  stimulated  much  activity  on  this  topic  in  the 
Operations  Research  community.  This  section  does  not 
comprehensively  review  all  works  that  have  been  written 
In  the  literature,  but  It  presents  a  brief  overview  of 
the  utilization  of  VRTs  in  simulation.  The  purpose  of 
this  section  Is  to  summarize  some  of  the  studies  of  the 
general  applicability  of  VRTs  in  simulation. 


Hammersley  and  Handscomb  (-1964)  reviewed  many  of 
the  simplest  ideas  of  variance  reduction  in  simple 
Monte  Carlo  problems  as  they  can  be  applied  in  the 
fields  of  Mathematics  and  Physics.  Their  most  easily 
understood  examples  and  outstanding  successes  were  the 
evaluations  of  integrals  and  applications  to  particle 
physics.  Handscomb  (1969,  p.252-262)  later  suggested 
that  VRTs  be  adapted  to  simulation.  He  acknowledged 
difficulties  in  predicting  the  effectiveness  of  the 
techniques  in  particular  situations,  but  he  did 
propose,  in  practice,  "...  to  proceed  by  more  or  less 
inspired  trial  and  error,  learning  by  experience  which 
tools  serve  one  best  [or  which  techniques  are 
effective]."  He  also  stated  that  it  may  be  much  harder 
to  tell  how  much  variance  reduction  may  occur  in  large 
and  complicated  simulation  problems.  These  issues 
remain  major  concerns  for  one  using  VRTs  in  large, 
complex,  stochastic  simulation  models. 

Moy  (1969,  pp .  263-288)  adapted  several  VRTs  to 
simulation  and  investigated  their  applicability  to 
queuing  systems.  He  concluded  that  VRTs  were  indeed 
capable  of  working  in  the  simulation  of  queuing 
systems.  Kleijnen  (1974,  Ch .  Ill),  who  has  written 
probably  the  most  comprehensive  and  most  referenced 
documentation  on  the  subject  of  VRTs  in  simulation, 
discussed  the  relevant  differences  between  sampling  in 
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Monte  Carlo  problems  and  sampling  in  stochastic 
simulation  models.  He  showed  that  VRTs  may  be  adapted 
to  accommodate  these  differences.  Kleijnen  also 
presented  a  detailed  description  and  critical  appraisal 
of  six  techniques  so  devised  for  utilization  in  the 
simulation  of  large  complex  systems.  These  VRTs  are 
stratified  sampling.  importance  sampling,  selective 
samp  ling,  control  variates,  antithetic  variates,  and 
common  random  numbers.  These  six  sampling  techniques 
have  become  the  most  well-known  and  popular  VRTs  in  the 
literature . 

Other  less-known  VRTs,  however,  have  been  applied 
to  simulation.  McGraft  and  Irving  (1974)  survey  some  18 
different  techniques  for  implementation  in  large  scale 
simulation  problems.  McGraft  and  Irving  include  a 
comprehensive  listing  of  the  characteristics, 
advantages  and  disadvantages,  and  criteria  for 
applicability  to  large  simulation  models,  and 
demonstrate  the  effectiveness  of  several  of  these 
techniques  with  a  military  simulation  application. 

Many  other  articles  and  papers  have  been  written  on 
the  subject  of  VRTs.  There  are  too  many  of  them  to  list 
in  this  thesis,  but  the  survey  ranges  from  specific 
techniques  to  more  general  methods  in  simulation 
experimentation.  Some  of  the  most  recent  papers  written 
about  the  general  applicability  of  VRTs  in  simulation 


are  Nelson  (1985)  and  Cheng  (1986).  Textbooks  that 
illustrate  the  applications  with  simple  but  excellent 
examples  of  variance  reduction  in  simulation  are  Gaver 
and  Thompson  (1973,  Ch .  12),  Fishman  (1978,  Ch .  3),  Law 
and  Kelton  (1982,  Ch.  11),  Morgan  (1984,  Ch .  7),  and 
Brat  ley,  Fox,  and  Schrage  (1987,  Ch .  2). 

C.  DESCRIPTION  OF  VRTs  USED  IN  THIS  THESIS 

Moy  (1969,  p.  263-288)  experimented  with  antithetic 
variates  and  stratified  sampling  and  showed  that  they 
are  indeed  capable  of  significantly  decreasing 
variability  in  the  simulation  of  simple  queuing 
problems.  Likewise,  we  wish  to  achieve  similar  results 
when  we  apply  them  to  the  simulation  of  our  stochastic 
combat  model.  We  do  this  in  Chapter  V.  In  this  section, 
we  discuss  the  underlying  conditions  and  fundamental 
concepts  in  the  applicability  of  antithetic  variates 
and  stratified  sampling  in  simulation. 

1 .  Antithetic  Variates 

The  method  of  antithetic  variates  is  relatively 
well-known  in  the  literature  of  variance  reduction  in 
simulation  (Kleijnen,  1974).  It  is  one  of  the  most 
useful  VRTs  because  of  its  simplicity  and  general 
applicability.  When  the  method  of  antithetic  variates 
is  used,  the  sampling  process  is  modified  by  the 
manipulation  of  random  numbers.  A  simulation  run 


produces  a  response  from  the.  -  original  sequence  of 
random  numbers  <  r  > ;  then,  a  second  simulation  run 
produces  an  antithetic  response  from  the  sequence  of 
the  complementary  random  numbers  < 1 -r  ) .  The  average  of 
the  two  responses  is  an  observation  on  the  sample 
output  data  of  the  stochastic  simulation  model.  The 
mean  value  of  this  sample  is  estimated  as  the  parameter 
of  interest. 

The  variance  of  this  estimate  is  reduced  if  the 

responses  of  the  first  and  antithetic  runs  of  each 

replication  are  negatively  correlated.  Besides  the 

interchanging  of  the  random  numbers  in  each  run,  two 

other  conditions  must  occur  to  produce  negative 

correlation  between  the  runs.  First,  each  response  must 

be  a  monotonic  function  of  its  respective  random  number 

stream;  that  is,  large  values  in  each  stream  of  random 

numbers  should  have  an  opposite  effect  on  the  response 

than  the  small  values,  and  vice  versa.  The  second 

condition  is  that  the  responses  to  the  events  in  the 

first  run  must  be  synchronized  with  the  responses  to 

the  events  in  the  antithetic  run.  Synchronization, 

defined  by  Kleijnen  (1974,  p.  193),  occurs 

...if  the  i ’ th  random  number  r^  generates  [in  the 
stream  of  the  first  run]  a  particular  event  (e.g., 
arrival  of  customer  d )  then  in  the  antithetic  run 
(1  -  r^)  should  generate  the  same  event  (i.e.,  not 

the  arrival  of  customer  d ’  where  d  ’  ^  d  and  not  a 
service  time). 
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If  the  antithetic  variates  methodology,  coupled 
with  required  conditions,  is  designed  into  the 
simulation,  then  the  average  of  the  two  negatively 
correlated  responses  will  tend  to  produce  an  estimate 
with  a  high  degree  of  precision.  That  is,  if  by  chance 
<r)  yields  a  response  above  the  value  of  the  true 
parameter  of  interest,  then  O-r)  should  yield  a 
response  below  the  value  of  the  true  parameter.  When 
these  responses  are  averaged,  the  deviations  between 
the  responses  and  the  true  parameter  approximately 
offset  each  other  resulting  in  relatively  small  net 
variability  in  the  output  data.  This  idea  can  be  shown 
mathematically.  Let  Xj  be  the  response  of  the  first 
run;  X2,  the  response  of  the  antithetic  run;  and  Y,  the 
average  of  X-j  and  X2. 

Y  -  (X,  +  X2)  /  2 

VARy  ■  1/4  *  {  VARxi  +  YARj(2  +  2  *  COV^i  ^2  ) 

=  1/2  *  (  1  +  C0RRxifX2^  *  V  ARx  i 

Clearly,  a  negative  C0RRxi,x2  reduces  VARy.  If 
C0RRxi(x2  equals,  or  is  close  to,  -1,  then  the  VARy  is 
mathematically  zero  or  very  close  to  it.  Hence,  the 
antithetic  sampling  is  designed  in  simulation  models  sc 
that  the  correlation  between  the  pair  of  responses  is 
as  close  to  -1  as  possible. 

1  9 


must  be 


.Monotonicity  and  synchronization 
designed  into  a  simulation  program  for  a  particular 
model  •  Kleijnen  (1974),  Law  and  Kelton  (1982),  and 
Bratley,  Fox,  and  Schrage  (1987)  are  excellent 
references  that  discuss  ways  to  do  this.  We  discuss  our 
design  to  achieve  these  two  conditions  for  antithetic 
variates  in  our  model  in  Chapter  IV.  As  stated  before, 
the  method  of  antithetic  variates  is  simple  to 
implement  and  requires  little  to  no  extra  computer 
time.  Because  of  simplicity  of  this  VRT ,  examples  of 
its  applications  are  illustrated  in  nearly  every 
textbook  that  considers  the  subject  of  VRTs. 

2  .  Stratified  Sampling 

The  stratified  sampling  technique,  discussed  in 
this  section  and  applied  to  the  simulation  model  in 
chapter  V  of  this  thesis,  is  a  different  version  of  the 
stratified  sampling  that  Moy,  Kleijnen  and  othe^ 
experts  on  VRTs  have  adapted  to  simulation.  Handscomb  ( 
1969,  p.  261)  calls  this  particular  version  of 
stratified  sampling  another  form  of  antithetic 
sampling.  Andrtasson  (  1972,  p.  6)  refers  to  it  as  an 
antithetic  transformation.  Gaver  and  Thompson  (1973, 
pp .  585-586)  name  it  stratification  extending  an 

antithetic  idea.  It  is  indeed  stratification  in  that 
the  sampling  process  is  modified  so  that  the  range  of 
random  numbers  is  divided  into  two  or  more  strata  from 
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which  the  simulation  runs  produce-  responses.  It  has  the 
antithetic  flavor  in  that  the  responses  in  all  strata 
are  averaged  together  to  get  an  observation  which  is 
part  of  the  samp]  -  '.\tput  data.  This  technique  is  also 
similar  to  antithetic  variates  in  that  its  estimator  is 
an  average  of  correlated  responses  (Gaver  and  Thompson 
1973,  p.  586).  Likewise,  this  estimator  tends  to  have  a 
smaller  variance. 

In  our  review  of  this  technique,  we  saw  no 
necessary  conditions,  like  those  for  the  antithetic 
variates,  for  this  technique  to  be  successful  in 
simulation.  The  design  of  stratified  sampling  into  our 
simulation  model  in  Chapter  IV  is  similar  to  t lat  one 
given  in  Gaver  and  Thompson  (1973,  p . 586 ) . 

D.  SUMMARY 

An  abundant  amount  of  material  has  been  written  on 
the  subject  of  variance  reduction.  Techniques  used  to 
reduce  the  variance  in  Monte  Carlo  problems  have  been 
adjusted  to  do  the  same  in  simulation  models.  The 
applications  of  VRTs  in  simulation  have  been 
illustrated  in  queuing  systems  and  simple  textbook 
problems  but  successful  applications  to  larger,  more 
complex,  real-world  stochastic  simulation  models  have 
not  been  so  amply  reported.  There  is  no  guarantee  that 
VRTs  will  work  spectacularly  for  every  situation  in  the 
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simulation,  and  when  they  do  work  it  is  necessary  to 
estimate  the  magnitude  of  the  variance  reduction.  Pilot 
tests  are  encouraged  to  help  resolve  these  issues. 
Antithetic  variates  and  modified  versions  of  stratified 
sampling  are  two  of  the  more  simple  and  easily  employed 
VRTs  and  will  be  applied  to  a  stochastic  combat  model 
in  Chapter  V. 
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III.  EFFICIENCY  OF  VARIANCE  REDUCTION 


A.  INTRODUCTION 

In  the  last  chapter  we  reviewed  some  studies  that 
involved  VRTs .  In  this  chapter  we  discuss  the  pr  b  era 
of  measuring  the  efficiency  of  a  VRT .  Comparing 
variances  of  a  parameter  of  interest  obtained  from  the 
simulations  with  and  without  the  use  of  a  VRT 
respectively,  on  an  ordinal  scale,  may  reveal  if  the 
VRT  works,  but  it  provides  little  information  about  how 
well  the  VRT  works.  Clearly,  a  quantitative  measure  is 
more  desirable.  Therefore,  the  manner  or  scale  on  which 
the  efficiency  of  a  VRT  is  measured  should  provide  as 
much  information  as  possible  on  the  performance  of  a 
VRT.  In  particular,  it  should  provide  at  least  some 
base  to  answering  the  following  questions: 

(i)  "Does  the  VRT  work?" 

(ii)  "If  so,  how  great  is  the  variance  reduction  in 
terms  of  increased  precision  for  estimating 
the  parameter  of  interest?" 

(iii)  "How  great  is  the  variance  reduction  in  terms 
of  simulation  time  saved  for  estimating  the 
parameter  of  interest?" 

(iv)  "What  are  the  tradeoffs,  if  any,  between  the 
potential  increase  in  precision  and  the  economy 
of  simulation  time  when  applying  VRTs?" 

In  the  next  section  we  examine  two  methods  that  are 
usually  used  in  the  literature  to  measure  the 
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efficiency  of  a  VRT .  We  evaluate  them  in  terms  of  how 
well  they  answer  the  questions  above.  In  the  third 
section,  we  offer  a  third  alternative  which  is  a  hybrid 
of  the  two  previous  methods  for  measuring  the 
efficiency  of  a  VRT.  This  third  method,  we  think, 
answers  all  four  questions  above  and  is  used  to  measure 
the  efficiencies  of  the  antithetic  variates  and  the 
stratified  sampling  techniques  whose  performance  is 
exhibited  in  this  thesis.  In  the  fourth  section  of  this 
chapter  we  show  how  to  use  the  third  method  of 
measuring  the  efficiency  of  a  VRT  to  obtain  the 
tradeoffs  between  increased  precision  and  reduced 
simulation.  The  last  section  is  a  summary  of  this 
chapter . 


B.  ASSESSMENT  OF  VARIANCE  REDUCTION 


In 

the  literature 

the 

efficiency  of  a 

VRT  is 

usual ly 

measured  by 

(1  ) 

a  decrease  in  the  variance 

( Method 

M 1  )  or  (  2  ) 

the 

relative  efficiency 

of  a 

simulation  to  obtain 

an 

estimate  using  a  VRT 

to  the 

efficiency  of  the  simulation  using  no  VRT  (Method  M2). 
Henceforth,  we  refer  to  a  simulation  without  the  use  of 
a  VRT  as  crude  simulation. 

Method  M 1  is  well  defined  in  Kleijnen  (1974,  pp . 
106-107).  Kleijnen  uses  this  method  by  defining  the 
efficiency  of  a  VRT  as  a  percentage  of  reduction  i-  variance 
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Method  #1  =  (Varg  -  Var^)  /  Var^)  *  1 00fc 


(2) 


where  Varg  and  Var^  are  variances  obtained  in  the  same 
amount  of  simulation  time  for  crude  simulation  and 
simulation  applying  a  VRT  respectively.  The  measure  of 
efficiency  of  a  VRT  which  Kleijnen  introduces  may  be 
interpreted  as  that  portion  of  the  variance  which  is 
not  achieved  by  crude  simulation  but  is  obtainable  in 
the  same  amount  of  simulation  using  a  VRT.  The  sign  of 
this  portion  determines  whether  the  VRT  increases  or 
decreases  the  precision;  a  positive  sign  reveals  an 
increase  and  a  negative  sign,  a  decrease.  The  magnitude 
of  the  portion  indicates  how  much  of  the  precision  is 
increased  or  decreased  respectively.  With  this  method 
we  can  also  see  that  the  VRT  has  an  identical  effect  on 
reducing  simulation  time  for  a  prescribed  level  of 
precision  as  it  does  on  increasing  precision.  Method  #1 
provides  answers  to  three  of  the  questions  stated  in 
the  lact  section,  but  it  does  not  resolve  the  question 
of  tradeoffs  for  increased  precision  and  reduced  time 
in  a  simulation  using  VRT. 

McGrath  and  Irving  (1974,  p.  295)  use  Method  #2  to 
measure  the  efficiency  of  a  VRT.  They  initially  used 
this  method,  shown  as  Equation  (3),  to  equate  the 


relative  advantage  gained  in  simulation  time  by  using  a 

\  i 

Method  M2  = 

Varg/Var^  *  (Simulation  Timegi )/( Simulation  Time^ )  (3) 

where  subscripts  0  and  1  represent  crude  simulation  and 
simulation  applying  a  VRT  respectively.  The  relative 
efficiency  that  McGrath  and  Irving  used  to  measure  the 
efficiency  of  a  VRT  results  in  a  factor  by  which  the 
efficiency  of  a  simulation  is  increased  or  decreased  by 
using  a  VRT.  If  the  value  of  this  factor  is  greater 
than  one,  then  the  VRT  works;  otnerwise,  it  does  not. 
The  magnitude  of  this  reduction  is  the  actual  value  cf 
the  factor.  For  example,  if  the  value  of  the  factor  is 
5,  then  the  simulation  applying  the  VRT  can  obtain  an 
estimate  in  1/5  the  simulation  time  re  ,uired  by  the 
crude  simulation  for  the  same  precision  level.  Method 
M2  may  be  viewed  either  as  the  reduction  in  simulation 
time  when  both  simulations  are  to  achieve  the  same 
variance,  or  as  an  increase  in  precision  when  both 
simulations  are  run  for  the  same  amount  of  time.  This 
method,  like  Method  M 1,  answers  only  the  first  three 
questions  proposed  in  the  first  section  of  this 


chapter . 


THE  HYBRID  METHOD 
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Methods  01  and  02  measure  increased  precision  at  a 
fixed  simulation  time,  or  a  reduced  simulation  time  at 
a  fixed  level  of  precision.  They  do  not,  on  the  other 
hand,  measure  increased  precision  at  a  level  of  reduced 
simulation  time,  or  vice  versa;  nor  do  they  provide  a 
means  .o  explore  such  a  possibility.  In  Chapter  I  we 
emphasize  that  variance  reduction  may  increase 
precision  and  reduce  simulation  time.  The  efficiency  of 
a  VRT ,  in  our  opinion,  should  reflect  both  effects  so 
that  we  can  explore  the  tradeoff  of  any  combination  of 
precision  and  simulation  time.  Method  03  offers  such 
possibility  and  answers  all  four  questions  in  the 
introduction  section  of  this  chapter.  It  is  a  mixture 
of  Methods  01  and  02  .  Method  03  has  Kleijnen’s  idea  of 
reduction  in  variance  and  McGrath  and  Irving’s  use  of 
relative  efficiency.  We  define  the  efficiency  of  a  VRT 
as  a  relative  efficiency  (RE),  as  shown  in  Equation  4, 
and  later  define  it  in  terms  of  increased  precision 
(IP)  and  reduced  time  ( RT ) . 


Method  03  = 

Efficiency-)  /  Efficiencyg  (4) 

where  Efficiency0  and  Efficiency-)  are  the  efficiencies 
of  th.  crude  simulation  and  simulation  applying  a  VRT 
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respectively.  The  definition  of  the  efficiency  of  a 
simulation,  identified  by  Equation  1,  is  the  inverse  of 
the  product  of  the  sampling  variance  of  the  parameter 
estimate  and  the  work.  Henceforth,  we  equate  work  to 
simulation  time,  which  is  the  total  time  of  the 
simulation  model  to  obtain  a  parameter  of  interest  and 
a  specified  variance.  Such  time  may  be  computed  as  n 
replications  times  T  (average)  simulation  time  per  run. 
If  k  runs  are  in  a  replication,  then  simulation  time 
equals  the  product  of  kn  runs  and  T  (average) 
simulation  time  per  run.  When  these  variables  are 
substituted  in  Equation  1 ,  the  efficiency  equation 
becomes  Equation  5: 

EFFICIENCY  -  1  /  (  Var  *  k  *  n  *  T )  (5) 

If  we  are  to  Increase  the  efficiency  of  a 
simulation  using  a  VRT ,  then  we  must  attempt  to 
decrease  one  of  the  parameters  in  Equation  2.  The 
variable  k  runs  per  replication  is  a  constant  of  the 
VRT.  Specifically,  the  antithetic  variates  constant  k 
is  two;  stratified  sampling  constant  k  is  three  in  our 
study  (it,  can  be  more);  and  for  no  VRT,  the  constant 
value  of  k  is  one.  The  variable  T  is  model  dependent; 
that  is,  its  value  depends  on  the  input  parameters  of 


the  model.  Attempts  to  decrease  this  variable  may  be 


futile;  futhermore  Bratley,  Fox,. and  Schrage  (1987,  p. 
48)  point  out  that  relatively  little  can  be  done  to 
decrease  it.  We  discussed  the  relationship  between  the 
variance  (Var)  and  replications  (n)  in  Chapter  I.  They 
are  ,in  essence, the  variables  we  wish  to  decrease. 
Throughout  this  thesis  we  interchange  the  phrases 
decrease  in  variance  with  increased  precision  and 
reduction  in  replications  with  reduced  time 
(simulation).  If  we  substitute  Equation  5  into  Equation 
4,  we  get  Equation  6.  Note  since  T  (average)  simulation 
time  per  run  is  the  same  for  both  simulations,  it  is 
left  out  of  the  equation. 

RE  =»  Varg/Vari  •  kg/^  *  n0/n-|  (6) 

If  the  RE  value  in  this  equation  is  greater  than 
one,  then  the  VRT  successfully  increases  the  efficiency 
of  the  simulation  and  is  said  to  be  strong;  otherwise, 
it  is  said  to  be  weak.  A  strong  VRT  decreases  the 
variance  so  that  precision  is  increased  and  simulation 
time  is  reduced.  A  weak  VRT,  on  the  other  hand,  does 
not  decrease  the  variance  as  well  as  a  strong  VRT;  in 
fact,  a  very  weak  (or  subversive)  VRT  may  increase  the 
variance,  which  causes  a  reauction  in  precision  and 
necessitates  an  increase  i,.  simulation  time.  In  most 
simulation  models  a  VRT  may  be  strong  for  certain 
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conditions  and  weak  for  other-  conditions.  In  this 
thesis,  we  look  for  such  characterizations  of  the 
antithetic  variates  and  the  stratified  sampling  in  the 
stochastic  combat  model  we  describe  in  the  next 
chapter.  In  the  next  section,  we  define  Equation(6)  in 
terms  of  increased  precision  and  reduced  simulation 
time . 


D.  TRADEOFFS  OF  GAINING  PRECISION  AND  SAVING  TIME 
Let  us  define  IP  as  a  decrease  in  variance, 

IP  -  ( Var0  -  Var x  )  /  Var0  (7) 

and  RT  as  a  reduction  in  simulation  time 


Then 


and 


RT  -  (n0k0  -  n^)  /  n0k0. 
Var^  -  (1.0  -  IP)  «  Var0 
nlkl  "  (1.0  -  RT )  x  n0k 


(8) 


(9) 

(10) 


If  we  substitute  Equations  (9)  and  (10)  into 
Equation  (6),  we  get  Equat ion ( 1 1 ) . 


RE  -  1/(1. 0  -  IP)  "  1/(1. 0  -  RT)  (11) 
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Equation  1 1  defines  the  relative  efficiency  which 
we  defines  as  Method  #3  of  measuring  the  efficiency  of 
a  VRT,  in  terms  of  increased  precision  and  reduced 
time.  This  equation  resolved  the  unanswered  question 
identified  as  (iv)  in  the  introduction  section  of  this 
chapter.  With  this  equation,  we  can  examine  any 
combination  cf  IP  and  RT.  For  example,  suppose  we 
measure  the  efficiency  of  a  VRT  to  have  a  RE  value  of  6 
for  the  same  amount  of  simulation  time  (Hint:  RT=0). 
Substituting  these  values  into  Equation  (11),  we  get  IP 
=  5/6  or  83.3%  increased  precision. 

Suppose  we  only  need  to  increase  the  precision  to 
75%  instead  of  83.3%,  then  we  can  substitute  the  values 
for  IP=3/4  and  RE=6  (RE  should  not  change)  into 
Equation  11.  We  now  get  RT»1/3  (Note,  we  increase  the 
precision  75%  and  reduce  the  simulation  time  33.3%). 
Likewise,  with  RE  =  6  for  the  efficiency  of  the  VRT, 
examples  of  other  combinations  are  ( IP*2/3 ,RT*1 /2 ) ; 

( IP* 1 /2 ,RT«2/3  )  ;  and  ( IP»0 ,RT»5/6 ) .  In  fact,  we  may  get 
any  combination  of  (IP,RT)  between  0  and  5/6.  Note, 
however,  if  we  want  to  increase  the  precision  beyond 
83.3%  or  5/6,  we  will  get  an  increase  in  simulation 
time.  That  is  the  tradeoff  in  terms  of  more  increased 
precision.  For  example,  we  will  have  to  increase  the 
simulation  time  to  2/3  or  66.7%  to  accommodate  an  IP  of 
90%  for  a  RE  value  of  6.  In  short,  the  information 
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obtained  from  method  03  is  that  ve  can  make  an  estimate 
more  precise  and  save  simulation  time  simultaneously. 

E.  SUMMARY 

In  the  literature,  there  are  generally  two  methods 
of  measuring  the  efficiency  of  a  VRT .  Method  01  is  a 
decrease  in  the  variance;  Method  02  is  the  relative 
efficiency  of  a  simulation  using  VRT  to  crude 
simulation.  Both  methods  may  determine  if  VRT  works  in 
a  simulation  model .  They  also  may  indicate  the 
magnitude  of  the  variance  reduction  in  terms  of  either 
increased  precision  for  a  fixed  simulation  time  or 
reduced  simulation  time  at  a  fixed  level  of  precision. 
In  this  chapter,  we  introduced  a  third  method  of 
measuring  the  efficiency  of  a  VRT.  It  is  a  hybrid 
between  Method  01  and  Method  02.  Method  03  offers 
exploration  into  the  tradeoffs  of  increasing  precision 
and  saving  simulation  time  for  any  efficiency  value  of 
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A.  INTRODUCTION 

In  the  last  chapter  we  discussed  how  to  measure  the 
efficiency  of  a  VRT  to  determine  how  much  we  may  save 
in  simulation  time  or/and  how  much  we  may  increase  the 
precision  of  a  parameter  obtained  by  crude  simulation. 
In  this  chapter  we  show  how  we  may  apply  VRTs  to  the 
simulation  of  a  combat  system.  The  combat  model  which 
we  have  chosen  to  simulate  ar.d  to  apply  the  VRTs  is  the 
BCD  Markovian  model  developed  in  the  doctoral 
dissertation  of  Abdul-Latif  Rashid  Al-Zayani  (1986).  A 
modified  version  of  this  model,  formulated  by 
Professor  Donald  P.  Gaver ,  is  in  Appendix  A.  This 
stochastic  model  may  seem  very  simple,  but  its 
simulation  provides  invaluable  insights  into  the 
applicability  of  VRTs  to  stochastic  combat  model. 

Beside  being  stochastic,  the  BCD  Markovian  model  is 
also  discrete  and  dynamic  in  nature;  hence,  it  is  a 
discrete-event  simulation  model.  We  refer  those  readers 
who  want  to  know  about  the  nature  of  discrete-event 
simulations  +  ^  simulation  textbook  such  as 
Morgan(  1  984  )  ,  L<  and  Kel ton(  1  982  )  ,  or  Bratley,  Fox, 


and  Schr age (  1 987  )  . 

In 

this  thesis,  we 

describe  the 

simulation  of  the 

BCD 

model  in  terms 

of  discrete 
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;  ;nts.  We  describe,  in  detail,  .the  crude  simulation  of 
the  BCD  model  in  the  next  section.  In  Section  C,  we 
show  how  we  applied  antithetic  variates  and  stratified 
sampling  to  this  simulation.  We  summarize  the  chapter 
in  the  last  section. 

B.  CRUDE  SIMULATION  OF  THE  BCD  MARKOVIAN  MODEL 

We  discuss  the  crude  simulation  of  the  BCD  model  in 
four  parts.  First,  we  describe  the  combat  scenario; 
second,  we  define  the  characteristics  of  the  model; 
third,  we  explain  the  simulation  of  the  combat  process 
in  the  model;  and  finally,  we  discuss  a  FORTRAN 
simulation  program  written  for  the  model. 

1 .  The  Combat  Scenario 

As  part  of  an  air  defense  command,  a  wing  of 
aircraft  defenders  is  responsible  for  defending  an  area 
against  a  hostile  air  attack  from  a  group  of  bombers. 
When  detection  of  an  incoming  threat  occurs  a  flight  of 
D  defenders  is  launched  to  engage  B  bombers  making  the 
attack.  When  the  two  groups  are  within  aerial  combat 
range  the  defenders  seek  a  one-to-one  combat  engagement 
with  the  bombers  at  a  rate  *.  Only  one  free  defender 
can  engage  a  free  bomber  in  combat;  a  bomber  will 
generally  attempt  to  avoid  any  engagement  with  a 
defender.  A  combat  engagement  lasts  until  either  the 
bomber  is  killed  or  the  defender  is  killed.  A  defender 
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kills  a  bomber  at  a  rate  a,  and  a  bomber  kills  a 
defender  at  a  rate  *.  Hence,  at  any  instant  during  the 
combat  process,  a  defender  is  either  free  and  searching 
to  fight  a  bomber,  fighting  a  bomber,  or  killed.  A 
bomber  is,  likewise,  either  free  and  eschewing 
engagement  with  a  defender,  engaging  in  combat  with  a 
defender,  or  killed.  The  combat  process  is  continued 
until  either  force  is  completely  killed  off  or  the 
duration  of  the  combat  period  is  terminated  after  T 
units  of  time. 

2 .  The  Characteristics  of  the  BCD  Model 

Hartman  (1985,  p.  2-18)  characterizes  the 
structure  of  a  combat  simulation  model  as  combat 
entitles ,  attributes ,  and  events .  We  use  these 
characteristics  to  simulate  the  combat  process  of  the 
BCD  mode1  in  the  next  subsection.  Combat  entities  in 
the  BCD  model  are  free  bombers,  free  defenders,  and 
combat  engagements .  Each  entity  has  attributes  that 
describe  a  combat  scenario.  For  the  bombers,  the 
attributes  are  the  number  of  bombers  and  the  rate  a 
bomber  shoots  down  a  defender;  for  the  defenders,  the 
number  of  defenders  and  the  rate  a  defender  shoots  down 
a  bomber;  and  for  the  combats,  the  number  of  combat 
engagements  and  the  rate  that  a  bomber  and  a  defender 
engage  in  combat. 
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Law  and  Kelton  (1982,  p . . 4 )  define  an  event  in 
a  discrete-event  simulation  as  an  occurrence  which 
changes  the  state  of  the  system.  The  BCD  model  has  five 
events.  The  first  event  is  the  initialization  of  the 
air  battle.  The  initialization  event  governs  the 
initial  battle  conditions.  The  next  three  events,  are 
the  interim  events  in  the  combat  process.  These  events 
are  (l)  a  combat  between  a  bomber  and  a  defender,  (2)  a 
defender  killing  a  bomber,  and  (3)  a  bomber  killing  a 
defender.  The  occurrence  of  an  interim  event  changes 
the  state  of  the  combat  process  at  time  t.  The  state  of 
the  combat  process  of  the  BCD  model  is  represented  by 
the  tri variate- Markov  process  { B( t ) , C( t ) , D( t ) ; t >0 ) ; 

where,  B(t)  is  the  number  of  free  bombers  at  time  t, 
C(t)  is  the  number  of  combat  engagements  at  time  t,  and 
D(t)  is  the  number  of  free  defenders  at  time  t.  As  a 
Markov  process,  the  combat  process  moves  from  state  to 
state  according  to  one-step  transit  probabilities  that 
depend  only  on  the  current  state.  The  fifth  and  last 
event  in  the  combat  process  is  the  termination  of  the 
air  battle.  The  termination  event  manifests  the  "end  of 
the  battle"  conditions.  The  values  of  the  state  of  the 
system  at  the  occurrence  of  the  termination  event 
reflects  the  battle  outcome.  These  values  are  the 
numbers  of  bombers  and  defenders  that  are  alive  at  the 
end  of  this  air  battle.  We  will  consider  these  values 
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as  the  parameters  of  interest  in- the  simulation  of  the 
BCD  model . 

3 .  Simulation  of  the  Combat  Process 

We  simulated  the  combat  process  by  maintaining 
a  "bookkeeping  account"  of  the  changes  in  the  state  of 
the  combat  process  as  the  events  occur .  The  process 
begins  in  the  initial  state  { B( t ) , C( t ) , D( t ) ; t  =  0  )  with 
the  initialization  event  being  the  commencement  of  the 
air  battle.  Henceforth,  we  let  a  value  of  B( t )  equal  b, 
a  value  of  C(t)  be  c,  and  of  D(t)  be  d.  The  interim 
events  change  the  value  of  the  state  of  the  combat 
process  as  following: 


EVENT 

New  combat 

Bomber  kills  Defender 
Defender  kills  Bomber 


STATE 


b- 1 , c+1 , d- 1 
b , c- 1 , d+1 
b+1 , c- 1 , d 


The  combat  process  spends  X(b,c,d)  units  of 
sojourn  time  in  state  (b,c,d)  until  another  event 
occurs.  The  sojourn  time  X(b,c,d)  is  a  random  variable 
distributed  exponentially  with  mean 

P(  b  ,  c  ,  d )  =  1  /  ( iBD  +  (a  +  #)C) 
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for  the  time  to  the  next  event  in- the  combat  process  of 


the 

air  battle, 

given 

that 

at  time 

t  the 

state  was 

(b,c 

,  v.  '  Equation 

12  is 

this 

soj  ourn 

time  . 

To  derive 

this 

expression , 

we  use 

the  inverse  transform 

method  to 

obtain  unit-mean  exponential  random  variables,  where  Vj 
is  the  jth  random  number  in  the  sequence  of  a  stream  of 
uniform  random  numbers.  The  inverse  transform  method  is 
discussed  in  the  simulation  textbooks  listed  in  the 
reference  section  of  this  thesis. 

X(b.c.d)  *  -  P( b , c , d )  *  ln( V j  )  (12) 

We  use  the  value  of  Equation  12  to  advance  the 
simulated  time  of  the  air  battle  as  indicated  by 
Equation  13. 

t  ■  t  +  X(b.c.d)  (13) 

The  combat  process  moves  to  another  state  when 
another  event  occurs.  The  probability  of  a  specific 
interim  event  occurring  is  governed  by  an  embedded 
Markov  chain  whose  transition  probabilities  are  as 
follows : 
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EVENT 


PROBABILITY 


New  Combat 

Defender  kills  Bomber 
Bomber  kills  Defender 


iBD  *  P(b,c,d) 
oC  *  P( b  ,  c , d  ) 
aC  *  P( b  ,  c , d ) 


We  again  use  the  inverse  transform  method  to 
obtain  the  conditions  for  an  interim  event  to  occur  and 
to  induce  the  change  in  the  state  of  the  combat 
process.  Vj  is  the  jth  random  number  in  the  sequence  of 
a  different  stream  of  random  numbers.  These  conditions, 
events,  and  changes  in  the  state  of  the  combat  process 
are  listed  below. 


CONDITION  EVENT 

V j  iiBD *  P(  b  ,  c  ,  d )  New  Combat 

V-.  >iBD» P( b  , c  ,  d )  Defender  kills  Bomber 
and 

V-j  <(iBD+<.C)*P(b,c,d) 

otherwise  Bomber  kills  Defender 


STATE 

b- 1 , c+1 ,  d+1 
b , c- 1 , d+i 

b  +  1  , c- 1  ,  d 


Thus,  we  (i)  generate  a  uniform  random  number 
to  choose  which  interim  event  has  occurred,  (ii)  update 
the  state  of  the  combat  process,  (iii)  generate  another 
uniform  random  number  and  transform  it  to  an 
exponential  random  variable  to  determine  the  unit  of 
time  until  the  occurrence  of  the  next  event  and  to 
advance  the  simulated  time  of  the  combat  process.  We 
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repeat  this  procedure  until  the  occurrence  of  the 
termination  event.  The  termination  event  occurs  when 
(1)  all  the  bombers  are  killed  (B(t)  +  C(t)  =  0)  or  all 
the  defenders  are  killed  (D(t)  +  C(t)  =  0)  or  (2)  the 

time  duration  of  the  air  battle  has  expired  (t  i T 
units  of  time).  At  the  end  of  the  aerial  battle,  the 
combat  process  is  in  state  (b,c,d).  from  which  we  can 
compute  the  number  of  live  bombers  B  (B(t)  +  C(t))  and 

the  number  of  live  defenders  D  (D(t)  +  C(t)).  These  are 
values  of  random  variables  for  one  possible  battle 
outcome . 


4 .  The  FORTRAN  Simulation  Program 

We  coded  the  crude  simulation  of  the  BCD  model 
in  FORTRAN.  This  FORTRAN  program,  consisting  of  a  main 
program  and  four  subroutines,  is  in  Appendix  B.  The 
main  program  begins  in  an  interactive  mode.  The  program 
reads  the  values  for  the  attributes  of  a  combat 
scenario  from  the  terminal  and  sends  them  to  the  BATTLE 
subroutine.  BATTLE  runs  N  replications  of  the  comb  vt 
process  and  returns  the  summary  statistics  of  t  i 
outcome  of  N  battles  to  the  main  program.  The  main 
program  sends  them  to  the  STAT  subroutine.  STAT 
analyzes  these  battle  statistics  in  terms  of  parameter 
estimates  and  returns  the  values  of  these  parameter 
estimates  to  the  main  program.  The  main  program  then 
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sends  these  values  of  parameter  estimates  to  a 
formatted  output  file. 

The  BATTLE  subroutine  calls  two  subroutines 
UNJ.FOR  and  EXPON  for  the  generation  of  U(0,1)  random 
numbers.  These  two  subroutines  implement  the 
congruential  pseudo-random  number  generator 


ui+1  *  168071^  mod  (2**31  -  1) 


(14) 


discussed  and  tested  by  Lewis  and  Orav  (1985,  Ch .  V). 
UNIFOR  generates  a  sequence  of  uniform  random  numbers 
for  the  selection  of  the  occurrence  of  an  interim 
event.  EXPON  generates  a  sequence  of  uniform  random 
numbers  for  the  computation  of  the  unit  of  sojourn  time 
in  a  state. 


The  STAT  subroutine  performs  statistical  output 
analysis  for  the  simulation.  It  computes  the  means  and 
variances  of  the  sample  distributions  of  live  bombers 
and  defenders.  The  sample  means  for  bombers  and 
defenders 


B  =  SUM  /  N 
5  *  SUM  /  N 


(15) 

(16) 


are  unbiased  (point)  estimators  of  E[B(t)]  and  E[D(t)] 
respectively.  Similarly,  the  variances 


4i 


B)«*2  /  N  *  (N-1  ) 
D ) *  *  2  /  N  x  (N-1  ) 


(17) 

(18) 


VARg  »  SUM  (  B±  - 
VARg  .  SUM  (  Dj.  - 

are  unbiased  estimators  of  VAR[E(B(t)]]  and 
VAR [E[D( t )] ]  respectively  (Larson,  1982). 

C.  PROGRAMMING  FOR  VARIANCE  REDUCTION 

In  Chapter  I,  we  noted  that  VRTs  modify  the 
sampling  of  random  numbers.  In  this  section,  we  discuss 
these  modifications  for  the  antithetic  variates  and 
stratified  sampling  in  the  simulation  of  the  combat 
process.  We  describe  the  changes  we  made  to  the  crude 
FORTRAN  simulation  program  for  the  simulations  using 
antithetic  variates  in  Section  1  and  stratified 
sampling  in  Section  2  respectively. 

1 .  Antithetic  Variates 

We  make  changes  to  the  subroutines  BATTLE, 
UNIFOR,  and  EXPON  of  the  crude  simulation  program  to 
use  the  antithetic  variates.  The  FORTRAN  program  for 
the  BCD  simulation  model  applying  antithetic  variates 
is  in  Appendix  C.  The  BATTLE  subroutine  computes  the 
values  of  the  parameters  for  one  replication  as  the 
average  values  of  the  battle  outcomes  from  a  pair  of 
runs  of  the  combat  process.  We  obtain  the  values  of  the 
battle  outcome  from  the  first  run  by  using  a  stream  of 
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uniform  random  numbers  <U)  and  those  of  the  second  run 
by  using  a  stream  of  complementary  random  numbers  (1- 
U).  Since  we  are  attempting  to  decrease  the  variance  of 
the  estimates  by  inducing  negative  correlation  between 
these  two  runs,  we  want  to  minimize  this  negative 
correlation.  We,  first,  induce  a  negative  correlation 
between  the  two  runs  by  creating  monotonici ty  between 
the  random  numbers  and  the  values  of  the  battle  outcome 
within  each  run.  We  then  minimize  this  negative 
correlation  by  synchronizing  the  sequences  of  random 
numbers  {  (J  )  and  the  complement  {  1 -U  >  (Brat  ley.  Fox, 
Schrage  1987,  p.47).  Kieijnen  (1974,  p.187)  shows  that 
a  random  variable  generated  by  the  inverse  transform 
approach  is  monotonic.  Hence  we  have  monotonicity  in 
the  simulation  since  we  used  the  inverse  transform 
method  to  generate  the  uniform  random  variables  in  the 
simulation  of  the  BCD  model. 

Law  and  Kelton  (1982,  p.  352)  indicate  that  the 
inverse  transform  approach  also  facilitates  the 
maintenance  of  synchronization  .  With  this  method,  we 
use  only  one  uniform  random  number  per  sequence  to 
obtain  the  desired  random  variable  for  each  event  in 
the  combat  process;  as  contrasted  with  other  methods, 
like  the  rejection  method,  where  we  may  use  many  random 
numbers  to  produce  a  single  value  for  the  desired 
random  variable  of  the  same  events.  Thus  we  initiate 
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synchronization  in  the  model  when  we  use  the  inverse 
transform  method;  nevertheless,  we  must  still  preserve 
it . 

We  risk  losing  this  synchronization  in  the  BCD 
model  because  the  number  of  interim  events  occurring  in 
the  combat  process  per  run  is  a  random  variable.  Hence 
the  number  of  events  occurring  in  the  combat  process  in 
the  first  run  may  not  be  the  same  as  the  number  of 
events  occurring  in  the  combat  process  in  the 
antithetic  run.  Consequently,  the  number  of  random 
numbers  needed  in  the  antithetic  run  generally  differs 
from  that  required  in  the  first  run.  This  phenomenon 
leads  to  the  random  number  (Uj)  in  the  first  run  not 
being  synchronized  with  the  random  number  ( 1 -U^  )  of  the 
antithetic  run  (Kleijnen  1974,  p.  193).  In  other 
words,  the  complement  of  the  J th  uniform  random  number 
< 1 -Uj  )  is  not  used  for  the  Jth  event  in  the  combat 
process  in  the  antithetic  run.  We  are  not  able  to 
control  the  random  number  of  interim  events  in  the 
combat  process,  but  we  can  manage  the  way  in  which 
UNIFOR  and  EXPON  generate  uniform  random  numbers  so 
that  synchronization  is  maintained  between  the  pair  of 
runs  per  replication. 

We  used  the  suggestions  of  Law  and  Kelton 
(1982,  p.352)  and  Bratley,  Fox,  and  Schrage  (1987, 
p.47)  to  maintain  the  synchronization  that  the  inverse 


transform  method  has  initiated  in  the  BCD  simulation 
model.  We  modify  subroutines  UNIFOR  AND  EXPON  to 
generate  separate  streams  of  random  numbers  <U)  and 
complements  (1-U)  before  simulating  a  pair  of  runs  of 
the  combat  process.  When  the  subroutine  BATTLE  calls 
UNIFOR  and  EXPON,  it  receives  from  each  a  two- 
dimensional  array  of  random  numbers,  where  the  first 
column  contains  the  stream  <U>  and  the  second  column 
contains  the  stream  <1-U).  Hence,  we  use  the  first 
column  for  the  first  run  and  the  second  column  for  the 
antithetic  run.  This  approach  guarantees  that  if  the 
j  th  event  in  the  first  run  uses  <Uj),  then  the  3 th 
event  in  the  antithetic  run  will  use  { 1 -Uj } .  We  do 
waste  some  of  the  random  numbers  in  the  arrays,  but  we 
do  it  Judiciously.  Since  the  number  of  random  numbers 
used  in  the  combat  process  is  a  random  variable,  we  use 
only  those  random  numbers  that  we  need  in  each  column 
and  throw  away  the  remaining  so  that  no  overlap  is 
possible  for  the  next  pair  of  runs.  As  a  result,  ve 
maintain  synchronization. 

The  last  change  we  make  to  the  crude  simulation 
for  the  utilization  of  the  antithetic  variates  is  in 
the  subroutine  BATTLE.  The  subroutine  BATTLE  computes 
the  values  of  parameters  for  each  replication  as 
f ol lows : 
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Bi  V 

(B1! 

+  bV)  / 

2 

(19) 

where 

B1  i  * 

the 

number  of 

live 

bombers 

from 

the 

first  run 

B2  i  * 

the 

number  of 

live 

bombers 

from 

the 

second  run 

and 

Dj.  *  ( D1  i  +  D2i)  /  2  (20) 

where 


*  the 

number 

of 

1  ive 

defenders 

from 

the  first  run 

*  the 

number 

of 

live 

defenders 

from 

the  second  run 

2 .  Stratified  Sampling 

As  we  stated  in  Chapter  II,  stratified  sampling 
resembles  the  antithetic  variates  procedures,  and  so  do 
the  changes  to  the  crude  simulation.  Hence  we  make 
changes  similar  to  those  in  the  simulation  using 
antithetic  variates  for  the  simulation  using  stratified 
sampling.  The  FORTRAN  program  for  the  BCD  simulation 
model  using  stratified  sampling  is  in  Appendix  D.  We 
modify  subroutines  UNIFOR  and  EXPON,  where  each 
generates  a  three-dimensional  array  of  uniform  random 
numbers  from  the  three  strata 


S,  -  (0,1 /3 )  , 


S2  -  ( 1 /3.2/3)  ,  S3  -  (2/3,1) 


before  simulating  three  runs  of  the  combat 
replication.  Note  this  does  not  need  to  be 
3;  we  could  have  done  more. 


process  per 
limited  to 
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Gaver  and  Thompson  (1973,  p . 586 )  describe  this 
approach.  First,  a  U(0,1)  random  number  u-|  1  is 
generated  and  placed  in  the  first  row  and  first  column 
of  the  three-dimensional  array.  Next,  1/3  is  added  to 
the  value  of  u^  ^  with  a  subtraction  of  one  if  needed  to 
get  U-J2  within  the  range  of  0  and  1.  u12  is  placed  in 
the  first  row  and  second  column  of  the  array.  Next,  1/3 
is  added  to  the  value  of  u12,  and  if  necessary 

subtracted  by  one,  to  get  u^.  u15  is  placed  in  the 

first  row  of  the  third  column  of  the  array.  If 

subroutine  BATTLE  calls  for  k  random  numbers,  then  k 
U ( 0 , 1  )  random  numbers  are  generated,  and  the  procedure 
obtains  a  value  for  each  of  the  kx3  cells.  BATTLE  uses 
the  first  column  of  random  numbers  in  the  array  for  the 
first  run,  the  second  column  for  the  second  run,  and 
the  third  column  for  the  third  run  of  the  combat 
process . 

The  values  of  the  parameters  for  each 

replication  are  the  average  values  of  the  battle 
outcomes  from  the  three  runs.  The  subroutine  BATTLE 
computes  the  values  of  these  parameters  as  follows: 

Bi  *  (B^  t  B2i  +  B^)  /  3  (21  ) 

where 

B1  =  the  number  of  live  bombers  from  the  first  run 
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B21 

B’l 

and 


where 


the  number  of  live  bombers  from  the  second  run 
the  number  of  live  bombers  from  the  third  run 


D1  * 

(D 

\  * 

D2i  +  D ^ 

)  /  3 

(22) 

the 

number 

of 

live 

defenders 

from 

the 

first 

run 

the 

number 

of 

live 

defenders 

from 

the 

second 

run 

the 

number 

of 

1  ive 

defenders 

from 

the 

third 

run 

D.  SUMMARY 

The  simulation  of  the  3CD  model  is  a  discrete-event 
simulation.  It  begins  with  the  initialization  event  and 
ends  with  termination  event.  The  simulation  of  the 
combat  process  involves  generating  a  sequence  of  U(0,1) 
random  numbers  to  select  interim  event  occurrences  with 
changes  in  the  state  of  the  process  and  generating 
another  sequence  of  U(0,1)  random  numbers  to  determine 
the  unit  of  time  until  the  next  event  occurs  and  to 
advance  the  simulated  time  of  the  combat  process. 

The  programming  of  the  antithetic  variates  and 
stratified  sampling  modifies  crude  simulation. 
Monotonicity  and  synchronization  are  required  in 
generating  the  uniform  numbers  for  the  simulation  using 
these  VRTs.  Generating  random  numbers  by  the  inverse 
transform  method  guarantees  monotonicity.  Generating 
sufficient  random  numbers  by  the  inverse  transform 
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method  and  in  multi-dimensional  arrays  initiates  and 
maintains  synchronization. 
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r .  APPLICATIONS  OF  THE  ANTITHETIC  VARIATES  AND 

STRATIFIED  SAMPLING 

A.  INTRODUCTION 

We  are  now  prepared  to  demonstrate  the  application 
of  the  variance  reducing  techniques  to  the  simulation 
of  a  combat  stochastic  model.  In  this  chapter  we 
illustrate  the  performance  of  the  antithetic  variates 
(AV)  and  stratified  sampling  ( SS )  in  the  simulation  of 
the  BCD  model.  In  Chapter  IV  we  stated  that  the  mean 
and  variance  of  the  parameters  of  interest  estimated 
from  simulation  are  used  to  analyze  the  output  data  of 
the  model.  Usually  the  estimated  mean  is  of  primary 
interest  to  decision  makers,  and  the  estimation  of  the 
variance  is  secondary.  Since  we  use  the  variance  of  the 
parameters  estimated  from  the  simulation  of  the  BCD 
model  to  exhibit  variance  reduction,  we  will, 
henceforth,  focus  on  the  variance. 

We  examine  the  applicability  of  AV  and  SS  in  the 
BCD  model  by  simulating  many  scenarios  of  the  air 
battle  and  recording  increases  in  simulation 
efficiency.  We  investigate  AV  and  SS  performance  by 
mapping  a  response  surface  that  characterizes  the 
efficiency  of  variance  reduction  in  the  model.  In  the 
next  section  we  specify  the  scenarios  and  discuss  the 
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application  of  using  AV  and  SS  in  the  simulation  of 
these  scenarios.  We  then  build  response  models  that 
describe  the  performance  of  the  AV  and  SS  for  these 
scenarios  and  discuss  their  experimental  results  in 
Section  C.  We  present  a  brief  summary  of  the  chapter  in 
the  final  section. 

B.  APPLICABILITY  OF  ANTITHETIC  VARIATES  AND  STRATIFIED 

SAMPLING  TO  THE  BCD  SIMULATION  MODEL 

In  Section  B  of  the  previous  chapter  we  described 
the  general  scenario  of  the  BCD  model.  In  this  section 
we  specify  various  combat  scenarios  to  observe  how  AV 
and  SS  are  applicable  to  the  simulation  of  the  BCD 
model.  Recall  that  seven  attributes  characterize  a  BCD 
scenario.  Since  a  change  of  the  values  of  one  of  these 
attribute  will  produce  a  different  scenario,  we  chose 
to  change  the  values  for  three  attributes.  We  simulated 
10,  30  and  50  defenders  against  10,  30,  and  50  bombers 
at  "end  of  battle"  times  of  25,  75,  and  125.  We 
maintain  the  a  rate  for  a  defender  killing  a  bomber  at 
.05,  the  e  rate  for  a  bomber  killing  a  defender  at  .05, 
and  the  t  rate  for  a  bomber  and  a  defender  entering  a 
combat  engagement  at  .005.  We  also  initialize  every 
simulation  with  zero  combat  engagements.  Thus,  we 
observe  the  responses  of  the  simulations  at  3  "end  of 
battle"  times  in  9  different  scenarios.  This 
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arrangement  comprises  a  set  of  27  independent 
simulations . 

We  run  this  arrangement  for  crude  simulation, 
simulation  using  AV ,  and  simulation  using  SS.  As  a 
result,  we  perform  a  total  of  81  different  simulation 
experiments.  In  order  to  make  a  fair  assessment  of  the 
applicability  of  AV  and  SS,  we  examine  the  variance 
obtained  from  the  same  amount  of  simulation  or  the  same 
numbers  of  simulated  battle  runs  for  every  simulation. 
We  run  90  battles:  this  equates  to  90  replications  in 
crude  simulation,  45  replications  in  simulation  using 
AV ,  and  30  replications  in  simulation  using  SS .  Table 
E.i  of  Appendix  E  contains  the  statistical  output  for 
crude  simulation;  similarly,  data  in  Tables  E.2  and  E.4 
are  from  simulations  using  AV  and  SS  respectively.  We 
use  Equations  6  and  7  to  measure  the  efficiency  of  the 
variance  reduction  (RE)  and  the  increase  in  precision 
of  the  parameters  from  the  simulation  (IP)  and  place 
the  AV  results  in  Table  E.3  and  SS  results  in  Table 
E.5. 

The  values  in  Tables  E.3  and  E.5  show  that  AV  and 
SS  respectively  are  applicable  in  the  BCD  simulation 
model.  A  RE  value  greater  than  unity  indicates  that  the 
VRT  is  effective  in  increasing  simulation  efficiency.  A 
positive  IP  value  exhibits  their  effectiveness  to 
increase  precision  of  the  desired  parameter.  With  these 
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two  values,  we  may  also  find  the  tradeoff  of  saving 
simulation  time  or  gaining  precision.  We  acknowledge 
that  these  values  are  all  values  of,  or  realizations 
of,  random  variables,  but  we  believe  that  the  tables 
show  that  the  variance  reduction  adheres  to  a 
stochastic  pattern.  That  is,  the  random  variables 
obtained  under  certain  scenarios  will  tend  to  have  the 
same  relationship  to  the  random  variables  obtained 
under  other  scenarios.  For  example,  the  data  in  the 
tables  suggest  that  high  RE  and  IP  values  correspond 
to  the  scenarios  that  start  combat  with  same  numbers  of 
bombers  and  defenders.  The  RE  and  IP  values  obtained 
under  these  scenarios  appear  consistently  higher  than 
the  values  under  all  other  scenarios.  Hence,  the 
variance  reduction  measured  by  these  RE  and  IP  values 
are  stochastically  greater  than  the  variance  reduction 
obtained  from  any  other  scenario.  Since  such  even 
combats  (i.e.  equal  combat  power)  are  inherently  more 
variable  in  outcome,  the  fact  that  variance  reduction 
is  greatest  there  is  certainly  welcome.  In  the  next 
section  we  attempt  to  conduct  a  more  thorough 
investigation  of  these  phenomena  so  that  we  may 
understand  how  the  AV  and  SS  perform  in  the  simulation 
of  the  BCD  model . 
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C.  PERFORMANCE  OF  ANTITHETIC  VARIATES  AND  STRATIFIED 
SAMPLING  IN  THE  BCD  SIMULATION  MODEL 

In  the  previous  section  we  saw  that  AV  and  SS  are 
applicable  to  the  simulation  of  the  BCD  model.  In  this 
section  we  examine  the  variability  of  uncertainty  in 
the  model  and  then  evaluate  the  applicability  of  AV  and 
SS  to  reduce  this  uncertainty.  We  explore  the  changes 
in  the  AV  and  SS  performance  and  examine  the 
relationships  of  factors  that  affect  these  changes. 
Results  of  this  analysis  reveal  the  characterization  of 
the  AV  and  SS  performance  in  the  BCD  model. 

1  .  Experimental  Design 

We  use  the  data  we  generated  in  Appendix  E  to 
fit  response  surfaces  that  describe  the  uncertainty  in 
the  values  of  parameters  in  the  BCD  model  and 
cKm:  acterize  the  performance  of  AV  and  SS  over  the 
prescribed  range  of  values  in  the  three  factors: 
initial  numbers  of  Bombers  and  Defenders  and  "end  of 
battle"  Time.  We  code  the  three  factors  as 

Xi  =  (Time  -  75)  /  50, 

Xg  =  (Defender  -  50)  /  20, 

X3  =  (Bomber  -  50)  /  20. 

Each  factor  has  5  levels.  Thus,  we  may  use  a  5: 
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factorial  design  to  fit  the  data  with  the  response 
surface  equation 

T3(y)=  00  +  »iXi  +  #2x2  +  #3X3  +  ®11x1**2 
+  92 2X2  *  * 2  +  e33x3  *  +  *  1  2 x  1  *x2 

+  o 1 ^ x i  *x^  +  923x2  * x3  (  20  ) 

where 

E(.y)  =  expected  response 
80=  intercept 

9^=  linear  coefficient  for  factor  i 
8 i i =  quadratic  coefficient  for  factor  i 
interaction  coefficient  for  the 
interaction  of  factors  i  and  j 
x^=  level  of  factor  i. 

We  seek  to  obtain  the  maximum  information  from 
every  observation;  therefore,  we  chose  a  3^  Fractional 
factorial  design.  This  design  is  the  cuboctahedron  plus 
three  center  points  (John,  1971).  The  three  center 
points  provide  an  unbiased  estimate  of  error  and 
repeated  observations  which  permit  us  to  test  for  Lack 
of  Fit  of  the  response  surface  equation  we  obtained.  We 
use  the  cuboctahedron  design  to  fit  data  for  three 
response  surfaces:  (1)  variability  of  uncertainty 
inherent  in  the  battle  outcomes,  (2)  efficiency  of  AV  , 
and  (3)  efficiency  of  SS .  This  design,  with  its  three 
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center  points,  and  the  data  to  which  it  is  used  to  fit 
are  shown  in  Tables  5.1,  5.2,  and  5.3.  The  variances 


data  in  Table  5.1  is  the  variability  of  uncertainty 
inherent  in  the  battle  outcome.  We  obtain  this  data 
from  the  variance  of  the  estimate  tabulated  in  the 
crude  simulation  table  in  Appendix  E.  The  RE  values 
data  in  Table  5.2  is  the  efficiency  of  AV  for  the 
estimation  of  the  defender  and  bomber  parameters.  We 
took  this  data  from  the  RE  values  in  Table  E.3.  In 
Table  5.3  is  the  RE  values  data  for  the  efficiency  of 
SS .  This  data  is  obtained  from  the  RE  values  in  Table 
E.5. 


VABLE  5.1  DESIGN  AND  VARIANCES  FOR  A  3  x  3  EXPERIMENT 
ON  THE  VARIABILITY  OF  UNCERTAINTY  INHERENT  IN  THE 
BATTLE  OUTCOME 


DECODED 

LEVELS 

CODED 

LEVELS 

VARIANCES 

Time 

Defender  Bomber 

x  1 

x2 

x3 

Defender 

Bomber 

25 

1  0 

30 

-1 

-1 

0 

.  0328 

.  0541 

125 

1  0 

30 

1 

-1 

0 

.  0036 

.  1  725 

25 

50 

30 

-1 

1 

0 

.  1687 

.0826 

1  25 

50 

30 

1 

1 

0 

.3339 

.0359 

25 

30 

1  0 

-1 

0 

-1 

.  0458 

.  0350 

1  25 

30 

1  0 

1 

0 

-1 

.  2247 

.  0045 

25 

30 

50 

- 1 

0 

1 

.  0628 

.1104 

125 

30 

50 

1 

0 

1 

.  0362 

.  3479 

75 

1  0 

1  0 

0 

-1 

-1 

.  0387 

.0358 

75 

50 

1  0 

0 

1 

-1 

.  1  427 

.  0063 

75 

1  0 

50 

0 

-1 

1 

.  0064 

.  1  654 

75 

50 

50 

0 

1 

1 

.  2717 

.2865 

75 

30 

30 

0 

0 

0 

,1618 

.  1895 

75 

30 

30 

0 

0 

0 

.  0925 

.0966 

75 

30 

30 

0 

0 

0 

.  1  693 

.  1  589 
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TABLE  5.2  DESIGN  AND  RE  VALUES  FOR  A3  x  3  EXPERIMENT 
ON  THE  EFFICIENCY  OF  THE  ANTITHETIC  VARIATES  IN  THE  BCD 
MODEL 


DECODED  LEVELS 

CODED  LEVELS 

RE  VALUES 

Time 

Defender 

•  Bomber 

x  1 

x2 

x3 

Defender 

Bomber 

25 

1  0 

’0 

- 1 

-1 

0 

2.5 

5  . 

4 

125 

1  0 

00 

1 

-1 

0 

1 .5 

3. 

9 

25 

50 

30 

-1 

1 

0 

6 . 0 

2  . 

2 

125 

50 

30 

1 

1 

0 

7 . 2 

1  . 

4 

25 

30 

1  0 

-1 

0 

-1 

3.9 

3  . 

,  0 

125 

30 

1  0 

1 

0 

- 1 

4 . 0 

1  , 

,  5 

25 

30 

50 

- 1 

0 

1 

2.8 

3  , 

.  6 

125 

30 

50 

1 

0 

1 

1  .8 

1  0 

.  9 

75 

1  0 

1  0 

0 

-1 

- 1 

8. 1 

6 

.  0 

75 

50 

1  0 

0 

1 

-1 

3.5 

1 

.  1 

75 

1  0 

50 

0 

-  1 

1 

1  .  1 

4 

.  2 

75 

50 

50 

0 

1 

1 

10.1 

9 

.  3 

75 

30 

30 

0 

0 

0 

6.4 

7 

.  0 

75 

30 

30 

0 

0 

0 

13.5 

7 

.  0 

75 

30 

30 

0 

0 

0 

10.8 

9 

.  2 

TABLE  5.3  DESIGN  AND 
ON  THE  EFFICIENCY 
MODEL 

RE  VALUES  FOR 
OF  STRATIFIED 

A3  X  3  EXPERIMENT 
SAMPLING  IN  THE  BCD 

DECODED 

LEVELS 

CODED 

LEVELS 

RE  VALUES 

Time 

Defender  Bomber 

x  1 

x2 

x3 

Defender 

Bomber 

25 

1  0 

30 

-1 

-1 

0 

2 . 6 

2 . 1 

1  25 

1  0 

30 

1 

-1 

0 

.  9 

2 . 2 

25 

50 

30 

- 1 

1 

0 

3.0 

1  .  9 

1  25 

50 

30 

1 

1 

0 

2 . 2 

1  .  2 

25 

30 

1  0 

- 1 

0 

- 1 

2 . 0 

2.8 

1  25 

30 

1  0 

1 

0 

-1 

2 . 6 

1  .  7 

25 

30 

50 

- 1 

0 

1 

2.9 

1  .  4 

1  25 

30 

50 

1 

0 

1 

1  .  7 

2 . 6 

75 

1  0 

1  0 

0 

-  1 

- 1 

1  .  6 

1  .5 

75 

50 

1  0 

0 

1 

- 1 

1  .  6 

1  .  2 

75 

1  0 

50 

0 

-1 

1 

2 . 0 

1  .  5 

75 

50 

50 

0 

1 

1 

3.0 

3.3 

75 

30 

30 

0 

0 

0 

3.7 

2.7 

75 

30 

30 

0 

0 

0 

2.9 

3 . 0 

75 

30 

30 

0 

0 

0 

2.5 

1  .5 
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2 .  Experimental  Analysis 

To  perform  the  statistical  analysis,  we  use  the 
Response  Surface  Regression  (RSREG)  procedure  in  the 
Statistical  Analysis  System  (SAS)  computer  software 
package  on  the  IBM  370  mainframe.  With  this  procedure 
we  were  able  to  obtain  a  second  order  response-surface 
equation  by  least-square  regression,  check  for  model 
adequacy,  test  for  lack-of-fit,  and  identify  critical 
surface  values  which  were  useful  in  helping  to  describe 
the  shape  of  the  surface. 

We  fitted  Equation  20  to  the  data  in  the 
respective  design  tables  in  the  previous  section  and 
obtained  multiple  response  surface  equations  and 
multiple  analysis  of  variance  (ANOVA)  tables  for 
corresponding  responses.  We  assess  the  adequacy  of  each 
equation  and  test  for  fit  from  its  corresponding  ANOVA 
table.  From  each  response  surface  equation,  we 
generated  additional  data  to  obtain  contour  plots.  We 
plotted  contours  of  variability  of  uncertainty  for  the 
initial  numbers  of  bombers  and  defenders  at  various 
Times.  Similarly,  we  plotted  contours  of  the 
efficiencies  of  AV  and  SS  for  initial  numbers  of 
bombers  and  defenders.  From  these  plot  we  were  able  to 
see  how  the  initial  numbers  of  bombers  and  defenders  in 
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combat  affect  the  variability,  of  uncertainty  in  the 
battle  outcomes  and  the  AV  and  SS  performance. 

a.  Variability  of  uncertainty  of  the  battle 
outcomes 

Equation  21  provides  an  adequate 
description  of  the  response  surface  that  characterizes 
the  variability  of  uncertain  in  the  expected  numbers  of 
live  defenders  at  the  end  of  the  battle  in  the  BCD 
model.  The  response  VDef>en(^er  is  the  expected  amount  of 
uncertainty  in  the  defender  estimate. 

.1412  +  .036xi  +  .104x2  —  *009x3  —  .014xi  *  *  2 
+  .049xiX2  +  .008X2**2  -  .051x1X3 
+  .040x2X3  -  .034X3**2  (21  ) 


TABLE  5.4  ANOVA  FOR  THE  EXPECTED  AMOUNT  OF  UNCERTAINTY 
IN  THE  DEFENDER  ESTIMATE 


SOURCE 

d.  f . 

SS 

MS 

F-RATIO 

Fitted  Surface 

9 

.  1  302 

.0145 

8 . 06 

Lack  of  Fit 

3 

.0102 

.  0034 

1  .  90 

Pure  Error 

2 

.  0036 

.  0018 

Total 

1  4 

.  1  440 

R- Square= .9043 

Mean  Var iance= . 1 1 94 

Std . 

Dev . = . 0525 

We  see  from  ANOVA  Table  5.4  that  the 
variation  in  the  variance  of  live  defenders  is 
insignificant  at  the  95%  level  (its  F-Ratio  of  8.06  is 
less  than  F# 9559,2  =  19.38).  The  Lack  of  fit  is  also 
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insignificant  (F-Ratio*  1.90  <  F.95;3,2  *  19.16). 

R-Square  value  is  .9043  which  indicates  that  about  90% 
of  the  variation  in  the  variance  of  the  expected  number 
of  live  defenders  is  accounted  for  by  Equation  21 . 
Further  analysis  reveals  that  the  response  surface  is 
shaped  like  a  rising  ridge.  The  plots  of  contours  at 
Figure  5.1  illustrate  the  nature  of  this  Response 
surface.  These  pictures  show  that  initial  numbers  of 
bombers  and  defenders  affect  the  variance  of  defenders 
at  Time  25.  At  Times  75  and  125  the  initial  numbers  of 
bombers  have  little  influence.  Here  the  variance  of  the 
expected  number  of  live  defenders  is  affected  solely  by 
the  increase  in  the  number  of  initial  defenders.  Hence, 
as  the  number  of  defenders  increases,  the  variability 
of  uncertainty  in  the  estimate  of  the  expected  number 

of  live  defenders  at  the  end  of  the  battle  increases. 

Equation  22  provides  an  adequate 
description  of  the  response  surface  that  characterizes 
the  variability  of  uncertainty  in  the  expected  numbers 
of  live  bombers  at  the  end  of  the  battle  in  the  BCD 
model.  V Bomber  is  the  expected  amount  of  uncertainty  in 
the  Bomber  estimate. 

vBomber=  -1483  +  .035xi  -  .002x2  +  •  1 04x^  -  .03lx-|**2 

-  .041XiX2  +  .032y.2##2  -  .067x-|X3 

+  .058x2X3  -  .007x3**2  (22) 
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Figure  5.1  Contour  Plots  of  the  Response  Surface  for 
the  Expected  Amount  of  Uncertainty  in  the  Defender 
Estimate 
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TABLE  5.5  ANOVA  FOR  THE  EXPECTED'  AMOUNT  OF  UNCERTAINTY 
IN  THE  BOMBER  ESTIMATE 


SOURCE 

d.f . 

SS 

MS 

F-RATIO 

Fitted  Surface 

9 

.  1  331 

.0149 

6.77 

Lack  of  Fit 

3 

.  0073 

.  0024 

1  .  08 

Pure  Error 

2 

.  0045 

.0022 

Total 

1  4 

.  1  408 

R-Squar e= . 9 1 88  Mean  Variance*. 1 ir3  Std.  Dev. =.0485 


ANOVA  Table  5.5  shows  that  the  variation  in 
the  variance  of  live  bombers  is  insignificant  at  the 
95%  level  (its  F-Ratio  of  6.77  is  less  than  F  = 

19.58).  The  Lack  of  fit  is  also  insignificant  (F-Ratio= 
1.08  <  F.95;3,2  =  19.16).  The  R-Square  value  indicates 
that  about  92%  of  the  variation  in  the  variance  of  live 
bombers  is  accounted  for  by  Equation  22.  Further 
analysis  reveals  that  this  response  surface  is  also 
shaped  like  a  rising  ridge.  The  plots  of  contours  at 
Figure  5-2  illustrate  the  nature  of  this  response 
surface.  These  figures  show  that  the  variance  of  the 
expected  number  of  live  bombers  generally  Increase  as 
the  initial  number  of  bombers  increases, 
b.  Antithetic  Variates. 

Equation  23  provides  an  adequate 
description  of  an  AV  response  surface  in  the  estimation 
of  the  expected  numbers  of  live  defenders. 
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VARIANCE  CF  BOMBER  ESTIMATE  AT  TIME  25 


^Defender®  10-20  -  .09xi  +  1  .70X2  -  *46x3  ~  4.25x1**2 

+  .55x^2  -  1.68x2*#2  -  .28X1X3 
+  3.40x2X3  -  2 . 85x3  *  *  2  (23) 

We  see  from  ANOVA  Table  5.6  that  the  variation  in  the 
efficiency  of  AV  is  insignificant  at  the  95%  level  (its 
F-Ratio  of  1.46  is  less  than  F .9559,2  *  19.38).  The 
Lack  of  fit  is  also  .'.nsignif  icant  (F-Ratio=.46  < 
F.95;3,2  =  19.16).  The  R-Square  value  is  .8497  which 

indicates  that  about  85%  of  the  variation  in  the 
Defender  RE  values  is  accounted  for  by  Equation  23. 
Here,  ^Defender  is  the  expected  simulation  efficiency 
of  AV  generated  to  reduce  the  uncertainty  in  the 
Defender  estimate. 


TABLE  5.6  ANOVA  FOR  THE  EXPECTED  EFFICIENCY  OF  AV  (RE 
VALUE)  GENERATED  TO  REDUCE  THE  UNCERTAINTY  IN  THE 
DEFENDER  ESTIMATE 


SOURCE 

d.f . 

SS 

MS 

F-RATIO 

Fitted  Surface 

9 

168.31 

18.70 

1  .46 

Lack  of  Fit 

3 

4.08 

1  .36 

.47 

Pure  Error 

2 

25  •  69 

12.84 

Total 

1  4 

198.08 

R-Square= .8497 

Mean 

RE  Value=5 . 54 

Std . 

Dev 

.=2.44 

Further  analysis  reveals 

that 

the 

response 

surface  is  shaped  like 

a  hill  with  a  gentle 

slope  on 

one  side  and  a 

fairly 

steep  slope 

on  the 

other  side. 

The  maximum  value  of  this  surface  occurs  in  the  BCD 


64 


scenario  that  begins  combat  with  50  defenders  and  40 
bombers  and  ends  the  battle  at  time  75.  The  plots  of 
contours  at  Figure  5.3  illustrate  the  nature  of  this  AV 
response  surface.  These  pictures  clearly  show  how  AV 
performs  in  different  scenarios  for  times  of  25,  75, 
and  125.  Beside  having  its  best  performance  in  a 
scenario  that  ends  at  time  75,  AV  appears  to  be  strong 
in  scenarios  that  initiate  the  air  battle  with  at  least 
30  defenders  and  30  bombers.  Its  weakest  performance 
seems  to  occur  in  those  scenarios  that  commence  combat 
with  no  more  than  30  defenders  and  30  bombers.  The 
plots  of  contours  show  that  the  efficiency  of  AV  is 
subversive  in  those  scenarios  whose  simulation 
initializes  the  air  battle  with  30  or  less  defenders 
and  40  or  more  bombers.  For  these  scenarios,  simulation 
efficiency  of  AV  may  often  increase,  instead  of 
decrease,  the  uncertainty  in  the  Defender  estimate.  We 
will  discuss  why  this  is  so  in  the  Experimental  Result 
Section.  Similar  analysis  of  the  AV  performance  is  made 
for  the  Bomber  RE  values.  An  adequate  description  of 
the  AV  response  surface  in  the  estimation  of  the 
bombers  is  characterized  by 

YBomber=  7-73  +  -44x1  "  -69x2  +  2-05x3  "  2.45xt**2 

+  .18x^2  -  2.05X2**2  +  2. 20X1X3 

+  2.50x2X3  -  . 53x3  *  *2  (  4  ) 
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EFFICIENCY  OF  AV  IN  DEFENDER  ESTIMATE  AT  TIME  25 


EFFICIENCY  OF  AV  IN  DEFENDER  ESTIMATE  AT  TIME  75 


EFFICIENCY  OF  AV  IN  DEFENDER  ESTIMATE  AT  TIME  125 


Figure  5.3  Contours  Plots  of  the  Responses  Surface  for 
the  Expected  Efficiency  of  AV  Generated  to  Reduce  the 
Uncertainty  in  the  Defender  Estimate 
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YBomber  is  the  expected  simulation  efficiency  of  AV  to 
reduce  the  uncertainty  in  the  Bomber  estimate.  The 
ANOVA  Table  5.7  indicates  that  the  proportion  of  the 
total  variation  in  the  Bomber  RE  values  accounted  for 
in  Equation  24  is  over  87%.  Furthermore,  this  variation 
is  insignificant  at  the  95%  level  (F-ratio  value  of 
8.19);  lack  of  fit  is  also  insignificant  (F-Ratio= 
1.61). 


TABLE  5.7  ANOVA  FOR  THE  EXPECTED  EFFICIENCY  OF  AV  (RE 
VALUE)  GENERATED  TO  REDUCE  THE  UNCERTAINTY  IN  THE 
BOMBER  ESTIMATE 


SOURCE 

d.f . 

SS 

MS 

F-RATIO 

Fitted  Surface 

9 

118.74 

13.19 

8.19 

Lack  of  Fit 

3 

14.17 

4.72 

1.61 

Pure  Error 

2 

3.23 

1  .61 

Total 

1  4 

136.14 

R-Square= .8722 

Mean  RE 

Value=5 . 05 

Std . 

Dev  .  =  1  .87 

Examining  this  response  surface  further  we 
find  that  the  shape  of  the  surface  changes  over  time. 
The  plots  of  contours  depicted  in  Figure  5.4  show  that 
the  shape  of  the  surface  looks  like  a  saddle  at  Time 
25,  a  gentle  slope  at  Time  75,  and  a  uniformly  rising 
ridge  at  Time  125.  The  critical  values  for  this  surface 
also  change  as  its  shape  changes.  Most  notable  are  the 
values  for  maximum  efficiency.  At  Time  25,  maximum 
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EFFICIENCY  OF  AV  IN  BOMBER  ESTIMATE  AT  TIME  25 


EFFICIENCY  OF  AV  IN  BOMBER  ESTIMATE  AT  TIME  75 


EFFICIENCY  OF  AV  IN  BOM3ER  ESTIMATE  AT  TIME  125 


Figure  b.4  Contour  Plots  of  the  Response  Surface  for 
the  Expected  Efficiency  of  AV  Generated  to  Reduce  tne 
Uncertainty  in  the  Eomber  Estimate 
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efficiency  occurs  in  scenarios  that  begin  fighting  with 
less  than  20  defenders  and  bomber.  By  Time  125,  maximum 
efficiency  has  shifted  to  the  scenarios  that  start 
with  at  least  40  defenders  and  bombers.  The  least 
amount  of  AV  reduction  occurs,  at  any  Time,  in  those 
scenarios  that  initialize  the  combat  simulation  with 
more  than  30  defenders  and  less  than  10  bombers. 

Here  is  a  summary  of  what  is  revealed  by 
the  above  analysis.  AV ,  in  general,  seems  to  be  the 
strongest  and  most  consistent,  and  equally-distributed 
between  closely-matched  pairs  of  bombers  and  defenders . 
Furthermore,  the  larger  the  evenly-matched  contest  the 
greater  the  variance  reduction.  When  the  defenders  and 
bombers  are  not  evenly  matched,  AV  is  not  as  consistent 
and  does  not  provide  equal  variance  reduction  in  the 
estimation  of  the  pair  of  parameters.  It  is  strong  in 
the  estimation  of  the  larger  combatants  and  weak  in  the 
estimation  of  the  smaller  ones. 

c.  Stratified  Sampling. 

We  analyze  the  efficiency  of  SS  in  the 
simulation  of  the  BCD  model  in  a  similar  manner  as  we 
analyzed  the  efficiency  of  AV .  If  we  analyze 

Equations  25  and  26  in  terms  of  ANOVA  Tables  5.8  and 
5.9  respectively,  we  will  get  results  similar  to  those 
we  obtained  in  the  last  section.  Therefore,  we  will 
forego  this  particular  analysis.  Note  that  Y^gf-gp^gp  is 
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the  expected  simulation  efficiency  of  SS  to  reduce  the 
uncertainty  in  the  Defender  estimate,  and  ^Bomber  is 
the  expected  simulation  efficiency  of  SS  to  decrease 
the  uncertainty  in  the  Bomber  estimate. 


YDef ender  = 

3.03  -  .39X! 

+  .49x2  + 

.  08x3  -  . 1 

5X-,  *  *2 

+  .  22x-)  x2  - 

7 0x2»*2  - 

.45x^3  + 

. 55x2x3 

-  . 58x3  *  *  2 

(25) 

TABLE  5.8 

ANOVA  FOR  THE 

EXPECTED 

EFFICIENCY 

OF  SS 

(RE 

VALUE)  GENERATED  TO  REDUCE  THE  UNCERTAINTY  IN  THE 
DEFENDER  ESTIMATE 


SOURCE 

d.  f . 

SS 

MS 

F-RATIO 

Fitted  Surface 

9 

8.24 

.92 

2.49 

Lack  of  Fit 

3 

.53 

.  18 

.47 

Pure  Error 

2 

.75 

.37 

Total 

14 

9-52 

R-Squar e= . 8666 1  Mean  RE  Value=2.27  Std.  Dev.  =  .51 


Y Bomber =  2-90  "  -06x1  ~  -03x2  +  .21x3  -  .41Xl*»2 

-  .20x^x2  -  .64x2**2  +  .  58x^3  +  .5 0x2x3 

-  . 36x3  *  *  2  (26) 
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TABLE  5.9  ANOVA  FOR  THE  EXPECTED  EFFICIENCY  OF  SS  (RE 
VALUE)  GENERATED  TO  REDUCE  THE  UNCERTAINTY  IN  THE 
BOMBER  ESTIMATE 


SOURCE  d.f.  SS  MS  F-RATIO 

Fitted  Surface  9  5.18  .58  19.33 

Lack  of  Fit  3  1.82  .61  20.33 

Pure  Error  2  .06  .03 

Total  14  7.06 

R-Square= . 7340  Mean  RE  Value=2.15  Std.  Dev.=.6l 


The  contour  plots  at  Figures  5.5  and  5.6  appear  to 
have  similar  features.  They  show  a  relatively  flat 
surface  except  at  the  corners.  The  corner  with  50 
Bombers  and  50  Defenders  has  the  highest  response  and 
the  other  corners  have  low  response.  These  plots 
suggest  that  the  SS  performance  is  generally  consistent 
in  all  but  a  few  scenarios  in  the  BCD  model.  Maximum 
efficiency  of  SS  occurs  in  those  scenarios  that 
initialize  simulation  with  equally  large  numbers  of 
bombers  and  defenders.  It  is  very  weak  in  those 
scenarios  that  begin  combat  with  either  less  than  10 
defenders  and  more  than  40  bombers  or  more  than  40 
defenders  and  less  than  10  bombers. 


3 .  Experimental  Results 

The  experimental  results  can  be  summarized  in 


Tables  5.10  and  5.11.  Table  5.10  shows  the 


relationships  between  the  AV  performance  and  the 
uncertainty  in  the  Defender  estimate  and  the  SS 
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EFFICIENCY  OF  SS  IN  DEFENDER  ESTIMATE  AT  TIME  25 


Figure  5.5  Contour  Plots  of  the  Response  Surface  for 
the  Expected  Efficiency  of  SS  Generated  to  Reduce  the 
Uncertainty  In  the  Defender  estimate 
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Figure  5.6  Contour  Plots  of  the  Response  Surface  for 
the  Expected  Efficiency  of  SS  Generated  to  Reduce  the 
Uncertainty  in  the  Bomber  Estimate 


performance  and  the  uncertainty  in  the  Defender 
estimate.  Similarly,  Table  5.11  shows  the  relationships 
between  the  AV  performance  and  the  uncertainty  in  the 
Bomber  estimate  and  the  SS  performance  and  the 
uncertainty  in  the  Bomber  estimate. 

We  further  examine  this  relationship  by 
analyzing  the  data  that  measure  the  uncertainty  (crude 
variance)  and  appropriate  variance  reduction  (RE 
Values)  in  Appendix  E.  After  we  applied  a  logarithmic 
transformation  to  the  data,  we  regress  the  RE  values 
on  the  crude  variance  data  and  observe  a  strong 
logarithmic  linear  relationship  between  uncertainty  and 
variance1  reduction. 


TABLE  5.10  RELATIONSHIP  BETWEEN  UNCERTAINTY  AND  THE 
EFFICIENCY  OF  VARIANCE  REDUCTION  IN  THE  DEFENDER 
ESTIMATE 


INITIAL 

UNCERTAINTY 

VARIANCE 

►.  DUCTION 

Def  enders 

Bombers 

Variance 

AV 

SS 

10 

10 

medium 

strong 

fair 

10 

30 

medium 

fair 

fair 

10 

50 

small 

weak 

weak 

30 

10 

large 

strong 

fair 

30 

30 

large 

strong 

strong 

30 

50 

medium 

fair 

fair 

50 

10 

1  arge 

strong 

weak 

50 

30 

large 

strong 

fair 

50 

50 

large 

s  t  rong 

fair 
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TABLE  5.11  RELATIONSHIP  BETWEEN  UNCERTAINTY  AND  THE 
EFFICIENCY  OF  VARIANCE  REDUCTION  IN  THE  BOMBER  ESTIMATE 


INITIAL 

UNCERTAINTY 

VARIANCE  REDUCTION 

Defenders  Bombers  Variance 

AV 

SS 

10 

10  medium 

strong 

fair 

10 

30 

1  arge 

strong 

fair 

10 

50 

large 

strong 

fair 

30 

10 

smal  1 

fair 

fair 

30 

30 

large 

strong 

fair 

30 

50 

large 

strong 

fair 

50 

10 

smal  1 

weak 

weak 

50 

30 

medium 

fair 

fair 

50 

50 

large 

strong 

fair 

This 

relationship  is 

manifested 

in  the 

mult iplicative 

equation  shown  in 

Table  5.12.  VAR 

1  s  the 

value 

of  uncertainty  or  variance  of  the  corresponding 

estimate  obtained  from  crude  simulation,  and  Y 

1  s  the 

simulation  efficiency  of  the  variance  reduction 

or  RE 

value 

for  the 

corresponding  estimate. 

TABLE 

5  .  12 

ANALYSIS  OF  THE 

RELATIONSHIP 

BETWEEN 

UNCERTAINTY  AND  EFFICIENCY  OF 

VARIANCE  REDUCTION  IN 

PARAMETER  ESTIMATES 

CORRELATION 

SET 

VRT 

ESTIMATE 

RELATIONSHIP 

COEFFICIENT 

POINT 

AV 

Defender 

Y-  10.43  *  VAR****. 

362  .8609 

. 00154 

AV 

Bomber 

Y-  9.73  "  VAR"". 

362  .8217 

. 00186 

ss 

Defender 

Y-  3.18  "  VAR"". 

144  .5601 

. 00032 

ss 

Bomber 

Y-  2.99  "  VAR"". 

147  .6054 

.  00058 
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The  correlation  coefficient  reveals  the 
strength  of  the  logarithmic  linear  relationships 
between  uncertainty  and  the  efficiencies  of  AV  and  SS . 
With  the  values  of  the  exponent  in  the  equations  being 
less  than  one  and  the  values  of  the  correlation 
coefficient  being  positive,  the  efficiencies  of  AV  and 
SS  are  observed  to  increase,  at  a  decreasing  rate,  as 
the  uncertainty  (variance)  increases.  We  obtained  the 
Set  Point  by  setting  Y=1  in  the  corresponding  equation 
and  solving  for  VAR.  At  this  value,  simulation 
efficiency  nether  increases  nor  decreases.  Now  if  we 
observe  a  value  of  uncertainty,  or  variance  obtained 
from  crude  simulation,  above  this  set  point,  then  we 
expected  to  get  an  efficiency  of  a  VRT  to  increase  the 
simulation  efficiency.  On  the  other  hand,  if  the  value 
is  below  the  set  point,  then  we  expect  the  efficiency 
of  the  VRT  to  decrease  the  simulation  efficiency. 

Here  is  the  bottom  line  on  AV  and  SS 

performance  in  the  BCD  model : 

1 .  If  we  apply  antithetic  variates  to  the  simulation 
of  the  BCD  model , 

a.  We  may  expect  the  variability  of  uncertainty 
in  the  defender  estimate 

(1)  to  decrease  if  the  variance  of  the 

estimate  obtained  from  crude  simulation 
is  at  least  .00154-,  and 
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(2)  to  decrease,  at  a  decreasing  rate,  with 
an  increase  in  the  variance  of  the 
estimate  obtained  from  the  crude 
simulation. 

b.  We  may  expect  the  variability  of  uncertainty 
in  the  bomber  estimate 

(  I  )  to  decrease  if  the  variance  of  the 
estimate  obtained  from  the  crude 
simulation  is  at  least  .00186,  and 

(2)  to  decrease,  at  a  decreasing  rate,  with 
an  increase  in  the  variance  of  the 
estimate  obtained  from  the  crude 
simulation. 

2.  If  we  apply  stratified  sampling  to  the  simulation 

of  the  BCD  model, 

a.  We  may  expect  the  variability  of  uncertainty 
in  the  defender  estimate 

( 1 )  to  decrease  if  the  variance  of  the 
estimate  obtained  from  crude  simulation 
is  at  least  .00032,  and 

(2)  to  decrease,  at  a  decreasing  rate,  with 
an  increase  In  the  variance  of  the 
estimate  obtained  from  the  crude 
simulation. 

b.  We  may  expect  the  variability  of  uncertainty 
in  the  bomber  estimate 

(1)  to  decrease  if  the  variance  of  the 
estimate  obtained  from  the  crude 
simulation  is  at  least  .00058,  and 

(2)  to  decrease,  at  a  decreasing  rate,  with 
an  increase  in  the  variance  of  the 
estimate  obtained  from  the  crude 
simulation . 

E.  SUMMARY 

We  illustrated  the  performance  of  the  antithetic 
variates  and  stratified  sampling  in  the  simulation  of 
the  BCD  model.  The  manifestation  of  their  pair 
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performance  was  characterized,  by  response  surface 
equations  and  plots  of  contour  lines.  Both  VRTs  were 
shown  to  be  effective  in  increasing  simulation 
efficiency,  but  they  perform  differently  in  the  BCD 
model.  AV  provides  the  largest  amount  of  variance 
reduction  but  is  more  volatile.  AV  increases  the 
simulation  efficiency  on  the  average  of  5  times  the 
crude  simulation;  it  is  strong  in  the  BCD  scenarios 
where  there  is  large  amount  of  uncertainty  in  the 
battle  outcomes  for  live  bombers  and  defenders,  and 
weak  in  those  scenarios  where  there  is  little  amount  of 
uncertainty  in  the  battle  outcomes.  SS ,  on  the  hand, 
has  a  more  consistent  performance.  SS  increases  the 
simulation  efficiency  at  a  mean  of  2  times  the  crude 
simulation.  It  performs  nearly  the  same  in  every 
scenario  except  where  the  uncertainty  is  large. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


A .  CONCLUSIONS 

The  objective  of  this  pilot  study  has  been  to 
investigate  the  applicability  of  VRTs  to  reduce  the 
inherent  variability  in  stochastic  combat  models.  We 
examined  the  effects  of  applying  AV  and  SS  to  the 
simulation  of  a  simple  stochastic  combat  model.  We  have 
now  shown  that  AV  and  SS  are  applicable  to  this 
stochastic  combat  model.  We  can  infer  that  these  VRTs 
are  Indeed  capable  of  working  in  stochastic  combat 
models,  and  their  prospects  in  larger  and  more  complex 
stochastic  combat  models  are  even  more  promising.  The 
conditions  of  monotonicity  and  synchronization  are 
essential  parts  of  the  design  of  the  simulation  program 
for  these  models.  Hence,  we  feel  that  sizable  Increase 
in  simulation  efficiency  is  possible  if  these 
requirements  are  met  in  the  simulation. 

The  experimental  results  of  applying  VRTs  to  the 
BCD  simulation  model  show  that  the  strength  of  AV  and 
SS  is  influenced  by  uncertainty.  A  strong  variance 
reduction  results  from  a  large  variance  of  the  estimate 
obtained  from  crude  simulation.  A  weak  variance 
reduction  is  caused  from  a  small  variance  of  the 
estimate  obtained  from  crude  simulation.  Hence,  sizable 
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and  consistent  variance  reduction  depends  on  large 
variability  of  the  simulated  output  from  the  stochastic 
combat  model.  Therefore,  the  variability  of  the  output 
from  larger  and  more  complex  stochastic  models  must 
also  be  large  enough  to  obtain  the  size  and  consistency 
of  simulation  efficiency  and  variance  reduction  one 
desires  from  the  applications  of  these  VRTs  to  such 
models . 

B.  RECOMMENDATIONS 

The  pilot  study  presented  in  this  thesis  provides 
a  base  for  further  studies  in  the  applications  of  AV 
and  SS  to  large-scale,  real  world,  stochastic  combat 
simulation  models.  Usually  complex  simulation  models 
have  many  subroutines  or  modules.  The  variability  of 
uncertainty  in  the  output  data  from  these  modules  may 
vary  from  low  to  high.  We  recommend  that  a  study  of 
this  matter  focus  on  the  degree  of  variability  of 
uncertainty  in  the  output  data  from  each  module.  The 
interest  of  the  study  should  be  concerned  with  the 
relationship  between  the  performance  of  the  VRT  and  the 
variability  of  the  output  data  from  each  module.  The 
results  should  indicate  whore  and  how  the  VRT  may  be 
used  in  the  model  in  order  to  maximize  the  simulation 
efficiency  of  the  model.  For  example,  if  the  study 
shows  that  a  VRT  perfo-ms  strongly  in  a  particular 
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module  whose  variability  of  output  data  is  large  and  it 
performs  poorly  in  another  module  whose  variability  of 
output  is  small,  then  the  study  should  recommend  that 
the  VRT  be  used  in  the  module  in  which  it  performs  best 
and  not  be  used  in  that  module  for  which  it  performs 
poorly.  Using  VRTs  in  the  module  with  small  variabi  lity 
would  most  likely  decrease  simulation  efficiency  for 
that  model  and,  at  worst,  suboptimize  the  overall 
performance  of  the  VRT  in  the  complex  combat  model . 
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APPENDIX  A 


FORMULATION  OF  THE  BCD  MARKOVIAN  MODEL 

by 

Pr  fessor  Donald  P.  Gaver 
No  al  Postgraduate  School 
Monterey,  California 

B  bombers  are  approaching  a  group  of  D  defenders. 
When  the  two  groups  approach  within  range  each  defender 
searches  for  a  bomber;  after  he  finds  one  they  engage 
in  combat.  Either  bomber  or  defender  may  win  the 
combat;  the  survivor  becomes  "free",  and  is  a  candidate 
for  the  next  combat.  In  general,  bombers  attempt  to 
avoic  combat,  defenders  seek  it  out. 

This  situation  becomes  a  tri-variate  Markov  chain 
if  the  following  state  is  defined:  { B( t )  ,  C ( t )  ,  D( t ) )  . 
Here  t  is  conveniently  measured  from  the  time  bombers 
and  defenders  are  close  enough  to  permit  combat  at  all, 
R(t)  is  the  number  of  free  bombers  at  +  lme  t 
thereafter;  ditto  for  D(t),  the  number  of  free 
defenders;  C ( t  )  Is  the  number  of  one-on-one  combats. 
Here  are  a  set  of  transition  rates: 

(1)  Combats  begin.  If  B(t)=b,  C(t)=c,  D(t)*d  and  if  * 
is  the  rate  at  which  free  bombers  are  found  by 
free  def^ders  then  with  probability 
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no  defender  finds  a  bomber  in  time  t.  Hence  the 
probability  that  a  defender  does  find  a  bomber  is 
ibdt  +  o ( t ) .  This  is  the  rate  at  which  free 
bombers  and  defenders  ge*  converted  to  combats: 
the  state 

(b,  c,  d)  «  (b-i  ,  c  + 1  ,  d-1  )  with  prob  ibdt. 

(2)  Defenders  win.  Same  initial  conditions.  If  a 
bomber  is  in  combat  with  a  defender  the 
probability  that  a  defender  shoots  down  the 
bomber  in  time  t  (combat  duration)  if  the  latter 
doesn’t  hit  the  defender  is 

P( Combat  duration  s  t  !  bomber  doesn’t  hit) 

=  1  -  e"at . 

Likewise,  the  probability  that  the  bomber  shoots 
down  the  defender  is 

P<Comtat  duration  i  t  !  defender  doesn’t  shoot) 

=  1  -  e_Bt. 
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Now  suppose  they  both  shoot »  doing  so 
independently.  Model  the  probability  that  both 
survive  to  t  as 

P<Combat  duration  >  t)  =  ( e_0M  (  )  * 

Now  the  probability  that  combat  lasts  until  t  and 
is  terminated  by  a  defender  shooting  down  the 
bomber  is 

P{ Combat  ends  is  ( dt ) ,  Defender  wins) 

=  ( e*atadt  )e“Bt 

=  <e-<«  +  •>t(a  +  #))(«/(«+•). 

This  shows  that  a  single  combat  duration  is 
exp(c  +  *)  and  the  event  of  a  defender’s  winning 
is  independently  o/(a  +  »).  Likewise,  the  combat 
duration  is  exp(o  +  »)  and  a  bomber’s  win  is, 
independently,  b/(  +  »).  If  there  are  c  combats 
going  on  then  the  first  combat  to  ends  does  so 
in  time  exp(c(a  +  »)). 

Hence 

(b,  c,  d)  ♦  (b,  c-1,  d-<-1  )  with  prob  ocdt 
(b,  c,  d)  ♦  (b+:,c-i,  d)  with  prob  cdt. 
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(3)  Simulation  and  Sojourn.  The  above  shows  that  we 
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may  simulate  the  combat  as  follows. 

(i)  You  are  in  state  (b,  c,  d).  Obtain  a 
sojourn  time  in  that  state  that  is 

exp(*bd  +  (a  +  b)c) 

i.e.  sbcd  =  1  /  (*bd  +  (a  +  o)c). 

The  system  stays  in  state  (b,c,d)  for 
time  [0,Sbcd). 

(ii)  With  probability 

ibd  /  (ibd  +  (a  +  #)c)  ♦  ( b  —  1  ,  c  +  1  ,  d-i  ) 

(NEW  COMBAT)  at  time  Sbcd. 

(iii)  With  probability 

ac  /  (ibd  +  (a  +  #)c)  ♦  (b,  c-1  ,  d+1  ) 

(DEFENDER  SHOOTS  DOWN  BOMBER)  at 
time  Sbcd. 

(iv)  With  probability 

«c  /  (»Ld  +  (o  +  t)c)  *  (b  +  1  ,  c- 1  ,d) 

(BOMBER  SHOOTS  DOWN  DEFENDER)  at 
time  Sbcd. 
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APPENDIX  B 

FORTRAN  PROGRAM  LISTING  FOR  THE  CRUDE  SIMULATION 

OF  THE  BCD  MODEL 


DIMENSION  BX( 1 01  ) ,DX(  1  01  ) ,SEED(2 ) , BOX(  2 ) ,DOX( 2 ) 

INTEGER  I.N.BB.DD.R 

REAL *4  X, Y,Z, BX.DX.TXT, BOX, DO X 

DOUBLE  PRECISION  SEED 

DATA  SEED  / 1 234 . 0D0 , 567890 1 23 . 0D0/ 

R=  2 
N  =  1  00 
C 

C  RECEIVE  INPUT  DATA  FROM  TERMINAL 
C 

WRITE( *  ,3  ) 

3  FORMAT ( 1 X . ’ENTER  THE  NUMBER  OF  BOMBERS’) 

READ  ’(12)’,  BB 

WRITE( * ,4  ) 

4  FORMAT( 1 X , ’ENTER  THE  RATE  WHICH  A  BOMBER  SHOOTS 
DOWN  A  DEFENDER ’ ) 

READ  ’ (F5 . 3 ) ’ ,  Y 
WRITE( * ,5 ) 

3  FORMAT ( 1 X , ’ENTER  THE  NUMBER  OF  DEFENDERS’) 

READ  ’ ( 12  )  ’  ,  DD 
WRITE( »  ,6) 

6  FORMAT ( 1 X , ’ENTER  THE  RATE  WHICH  A  DEFENDER’ 

&  ’SHOOTS  DOWN  A  BOMBER’) 

READ  ’ (F5.3)  ’  ,  X 

WRITE( * ,7  ) 'ENTER  THE  RATE  WHICH  FREE  DEFENDERS’ 

&  ’FIND  FREE  BOMBERS’ 

7  FORMAT ( 1 X , A ) 

READ  ’ (F5.3) ’ ,  Z 

WRITE( * ,8) ’ENTER  THE  TIME  DURATION  OF  THE 
&  ’BATTLE’ 

8  FORMAT ( 1 X , A ) 

READ  ’ (F5.3) ’ ,  TXT 

C 

C  RUN  REPLICATIONS  OF  N  BATTLES  AND  OBTAIN  SUMMARY  OF 
C  N  BATTLES 

C 

CALL  BATTLE( N , R , SEED , X , Y , Z , BB , DD , TXT , BX , DX ) 

C 

C  COMPUTE  STATISTICAL  OUTPUT  ANALYSIS  OBTAIN  PARAMETER 
c  ESTIMATES 
C 

CALL  ST AT (N,R,BX,DX,BOX,DOX) 

C 
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C  FORMAT  AND  PRINT  OUTPUT  OF  PARAMETER  ESTIMATES 
C 

WRITE(3,279)  N 

279  FORMAT( 15X , ’SAMPLE  SIZE '  ,  16  , / 15X  ,  6(  ’ - '  )  ,  IX  , 

&  4  (  ’  -  ’  )  ) 

WR ITE( 3 , 280  )  BB.DD 

280  FORMAT (//1X,I4,2X, ’BOMBERS’ ,6X, ’VERSUS’ ,6X, 
&  14 ,2X  ,  ’DEFENDERS ’  , / IX , 1 3 (  ’ - ’  )  ,  18X  ,  14 (  ’  -  *  )  ) 

290  WRITE  (  3,300  )TXT 

300  FORMAT ( / / 18X ,  ’TIME'  , F6 . 1  .  / . 4 7 (  ’ - ’  )  ) 

WR ITE( 3,310  ) 

310  FORMAT ( 1 9X , ’BOMBER’ , 2X , ’DEFENDER’ ,/18X, 

&  7(  )  ,2X,8(  ’-’  ) ) 

WRITE( 3,320)  BOX ( 1  )  , DOX ( 1  ) 

320  FORMAT { IX , ’AVERAGE’ , 7X.2F10.4 ) 

WR ITE( 3 , 330  )  BOX ( 2 ) , DOX ( 2  ) 

330  FORMAT ( IX, ’VARIANCE’ .6X.2F10.4 ) 

STOP 

END 

C 

C 

C  SUBROUTINE  BATTLE 
C 

SUBROUTINE 

&  BATTLE ( N , R , SEED ,X,Y,Z. BB.DD, TXT , BX . DX ) 
INTEGER  BB , DD , I ,K,N,R 

REAL  X ( N+l )  ,DX(N+1 ) , GA , SOJ , X , Y . Z , B . C . D . NC , 
&  BW.DW, INF, T.TIME.TXT 

DOUBLE  PRECISION  SEED(R) 

C 

C  INITIALIZE  STATISTICAL  COUNTERS 
C 

BX( N+l )-  0 . 0 
DX( N+l  )-  0 . 0 

C 

C  RUN  N  REPLICATIONS 
C 

DO  200  I-l.N 
C 

C  INITIALIZE  START-TO-BATTLE  VALUES 
C 

T-0. 0 
C-0. 0 

B-REAL( BB ) 

D-REAL( DD ) 

BW-  0.0 
DW-  0.0 
NC-  1.0 
C 

C  OBSERVE  OCCURRENCE  OF  AN  EVENT 
C 
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100  CALL  UNIFOR( SEED( I ) ,GA, 1 ) 

C 

C  DETERMINE  NEXT  INTERIM  EVENT  AND  UPDATE  THE  STATE 
C  (B.C.D) 

C 

IF( GA  .LE.  BW'  THEN 
B=B  +  1.3 
C-C  -  1.0 
D  =  D 

ELSE  IF  (CA  .LE.  ( BW+NC )  )  THEN 
B=B  -  1.0 
C  =  C  +  1.0 
D=D  -1.0 
ELSE 
B  =  B 

C=C  -  1.0 
D*D  +  1.0 
END  IF 

C 

C  COMPUTE  MEAN  TIME  IN  STATE  (B.C.D) 

C 

IF  ((B  .EQ.  0.0  .OR.  D  .EQ.  0.0)  .AND.  C 
&  .EQ.  0.0)  THEN 

INF=  1000000.0 

ELSE 

INF*  1  . 0/ ( Z*D*B  +  (X  +  Y ) *  C  ) 

END  IF 

C 

C  GENERATE  SOJOURN  TIME  IN  STATE  (B.C.D) 

C 

CALL  EXPON v 1  EED( 2  )  , SOJ , 1  ) 

TIME*  -INF  ALOG(SOJ) 

C 

C  ADVANCE  THE  SIMULATED  TIME  OF  THE  AIR  BATTLE 
C 

T-  T  +  TIME 
C 

C  COMPUTE  PROBABILITY  OF  NEXT  INTERIM  EVENTS  OCCURRING 
C 

BW-  Y  * C  *  INF 
DW-  X  *C  *  INF 
NC*  Z»B*D*INF 
C 

C  CHECX  CONDITIONS  FOR  OCCURRENCE  OF  TERMINATION  EVENT 


C 

C 

C 

C 


IF  (T  .  L'T  .  TXT)  GOTO  100 
ACCUMULATE  SUMMARY  OF  N  BATTLE  OUTCOMES 
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BX( I  )  =  B  +  C 
DX(  I  )■  D  +  C 
BX( N  +  1  )*  BX ( N+ 1  )  +  BX (  I ) 

DX ( N+ 1  )=  DX (N+1  )  +  DX (  I ) 

200  CONTINUE 
RETURN 
END 
C 
C 

C  SUBROUTINE  EXPON 
C 

SUBROUTINE  EX  PON ( SEED2 , A2 , K ) 

INTEGER  I.K 

REAL  STIM, A2 , INF 

DOUBLE  PRECISION  EFF , SEED2 

EFF=  2147485647. 0D0 

SEED2=DMOD( 1 6807 . 0D0  *  SEED2 , EFF ) 

A2=  SEED2/EFF 
RETURN 
END 
C 
C 

C  SUBROUTINE  UNIFOR 
C 

SUBROUTINE  UNIFOR( SEED1 , A1 , K ) 

INTEGER  I.K 
REAL  A 1 

DOUBLE  PRECISION  EFF.SEED1 

EFF=  21 47483647. 0D0 

SEED  1  =*DMOD(  1  6807. 0D0  *  SEED1.EFF) 

A 1 =  SEED1/EFF 

RETURN 

END 

C 

C 

C  SUBROUTINE  STAT 
C 

SUBROUTINE  STAT( N . R . BX . DX . BOX . DOX ) 

INTEGER  J.R.N 

REAL  BX ( N+1  ) ,DX(N+1  )  ,BOX(R) ,DOX(R) 

BOX ( 2  )  »  0.0 
DOX ( 2  )  *  0.0 

C 

C  COMPUTE  THE  ESTIMATES  OF  THE  SAMPLE  MEAN  AND 
C  VARIANCE 

C 

BOX (  1  )  *  BX( N  +  1  )/N 
DOX ( 1 )  =  DX( N+1 )/N 
DO  260  1= 1  ,N 

BOX ( 2  )  =  BOX ( 2  )  +  ( BX ( I ) -  BOX (1  )  )  *  *  2 
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DOX( 2  )  *  DOX (  2  )  +  ( DX( I ) -DOX ( 1  )  )  *  *2 
260  CONTINUE 

BOX( 2  )  *  BOX( 2)/(N*(N-1  )) 

DOX( 2  )*  DOX (2)/(N*(N-1  )) 

RETURN 

END 
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APPENDIX  C 

FORTRAN  PROGRAM  LISTING  OF  THE  BCD  SIMULATION 
USING  ANTITHETIC  VARIATES 


DIMENSION  BX( 1 01  ) ,DX( 1 01  ) , SEED( 2) , BOX( 2 ) ,D0X(2 ) 

INTEGER  I.N.BB.DD.R 

REAL  *  4  X ,Y ,Z , BX ,DX ,TXT .BOX ,DOX 

DOUBLE  PRECISION  SEED 

DATA  SEED  / 1 234 . 0D0 , 567890 1 23 . 0D0/ 

R  =  2 
N  =  50 

RECEIVE  INPUT  DATA  FROM  TERMINAL 
WR ITE( *  ,3) 

FORMAT ( 1 X , 'ENTER  THE  NUMBER  OF  BOMBERS’) 

READ  ’  ( 12  )  ’  ,  BB 
WRITE( *  ,4) 

FORMAT( 1 X , 'ENTER  THE  RATE  WHICH  A  BOMBER  SHOOTS’ 
&  ’DOWN  A  DEFENDER' ) 

READ  ’ ( F5 . 3  )  ’  Y 
WRITE( * .5) 

FORMAT ( 1 X ,  ENTER  THE  NUMBER  OF  DEFENDERS’) 

READ  ’  ( 12)  ’  ,  DD 
WRITE( * ,6 ) 

FORMAT( IX , ’ENTER  THE  RATE  WHICH  A  DEFENDER 
&  SHOOTS  DOWN  A  BOMBER’) 

READ  ’ (F3 .3) ’ ,  X 

WRITE( * ,7)  ’ENTER  THE  RATE  WHICH  FREE  DEFENDERS’ 
&  ’FIND  FREE  BOMBERS’ 

FORMAT ( 1 X , A) 

READ  ’ ( F5 . 3 ) ' ,  Z 

WRITE( * ,8) ’ENTER  THE  TIME  DURATION  OF  THE’ 

&  ’BATTLE’ 

FORMAT ( 1 X , A) 

READ  ’ (F5.3) ’ .  TXT 

RUN  REPLICATIONS  OF  N  BATTLES  AND  OBTAIN  SUMMARY  OF 
N  BATTLES 

CALL  BATTLE( N , R , SEED ,  X , Y , Z , BB , DD , TXT , BX , DX ) 

COMPUTE  STATISTICAL  OUTPUT  ANALYSIS  OBTAIN  PARAMETER 
ESTIMATES 

CALL  STAT( N , R , BX , DX , BOX , DOX ) 
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C  FORMAT  AND  PRINT  OUTPUT  OF  PARAMETER  ESTIMATES 
C 

WR ITE{ 3,279)  N 

279  FORMAT ( 1 5X , 'SAMPLE  SIZE ' , 16 , / 1 5X , 6( ' - ’ ) , 

&  1X,4(  ’-'  )  ) 

WRITE( 3 , 280  )  BB,DD 

280  FORMAT (//IX, 14, 2X, 'BOMBERS' ,  6X, 'VERSUS' , 

&  6X.I4.2X. 'DEFENDERS'  , / 1 X . 1 3(  ’ - ’  )  , 1 8X . 1 4(  ' - ’  )  ) 
290  WRITE  (3,300  )TXT 

300  FORMAT ( //18X ,  'TIME'  , F6 . 1  , /  ,  47 (  ’ - ’  )  ) 

WRITE( 3,310) 

310  FORMAT ( 1 9X ,  'BOMBER'  ,2X,  'DEFENDER'  , 

&  / 1 8X , 7 (  ’  =  ’  )  ,  2X , 8 (  *  =  ’  )) 

WRITE( 3,320)  BOX ( 1  }  , DOX ( 1  ) 

320  FORMAT ( 1 X , 'AVERAGE ’ ,  7X , 2F1 0 . 4 ) 

WRITE(  3,330)  BOX  (  2  )  ,  DO  .  t,  ) 

330  FORMAT( 1 X ,  'VARIANCE 6.  -FI  0.4) 

STOP 

END 


C  SUBROUTINE  BATTLE 
C 

SUBROUTINE 

&  BATTLE( N , R , SEED , X , Y  »  Z , BB . DD , TXT , BX , DX ) 
INTEGER  BB.DD.H, I , J .W.K.N.R 
REAL  GA( 2, 1  000) ,S0J(2, 1 000  )  ,BX(N+1  )  , 

&  DX ( N+1 ) , BAT( 50,2), DAT( 50,2), 

&  X, Y, Z, T, TXT, TIME, INF, NC.BW.DW.B.C.D 
DOUBLE  PRECISION  SEED(R) 

K=  1000 
C 

C  INITIALIZE  STATISTICAL  COUNTERS 
C 

BX( N  +  1  )=  0.0 
DX( N+1 )=  0.0 

C 

C  RUN  N  REPLICATIONS 

C 

DO  200  1=1 ,N 

CALL  SOJ OUR( SEED( 1  )  , SOJ , R , K ) 

CALL  STATE ( SEED( 2 ) , GA , R , K ) 

DO  175  J=1 , K 

C 

C  INITIALIZE  START-TO-BATTLE  VALUES 


H  =  0 
T  =  0 . 0 
C  =  0 .0 

B=REAL( BB ) 
D=REAL( DD ) 
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BW=  0.0 
DW=  0.0 
NC=  1 . 0 

OBSERVE  OCCURRENCE  OF  NEXT  EVENT 
00  H=H+1 

DETERMINE  NEXT  INTERIM  EVENT  AND  UPDATE  THE  STATE 
(B.C.D) 

IF( GA( J , H )  .LE.  BW)  THEN 
B=B  +  1.0 
C=C  -  1.0 
D  =  D 

ELSE  IF  ( GA( J , H  )  .LE.  ( BW+NC ) )  THEN 


B  =  B  - 

1  .  0 

C  =  C  + 

1  .  0 

D  =  D  - 

1 . 0 

ELSE 

B  =  B 

c=c  - 

1  .0 

D  =  D  + 
END  IF 

1  .  0 

COMPUTE  MEAN  TIME  IN  STATE  (B.C.D) 

IF  ((B  .EQ.  0.0  .OR.  D  .EQ.  0.0)  .AND.  C 
&  .EQ.  0.0)  THEN 

INF=  1000000.0 

ELSE 

INF  =  1 . 0/( Z*D*B  +  (X  +  Y)*C) 

END  IF 

COMPUTE  SOJOURN  TIME  IN  STATE  (B.C.D) 

TIME  =  -INF  *  ALOG( SO J ( J , H )  ) 

ADVANCE  THE  SIMULATED  TIME  OF  THE  AIR  BATTLE 
T=  T  +  TIME 

COMPUTE  PROBABILITY  OF  NEXT  INTERIM  EVENTS  OCCURRING 


BW=  Y  *C* INF 
DW=  X  *C  *  INF 
NC=  Z  *  B *D  *  INF 

CHECK  FOR  OCCURRENCE  OF  TERMINATION  EVENT 
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IF  (T  .LT.  TXT)  GOTO  100 
C 

C  RECORD  RESULTS  OF  BATTLE 

C 

BAT(  I ,  J  )  =  B  +  C 
DAT( I ,  J  )  =  D  +  C 
175  CONTINUE 

C 

C  ACCUMULATE  SUMMARY  OF  N  BATTLE  OUTCOMES 
C 

BX( I ) = ( BAT( 1,1)  +  BAT(I,2))*.5 
DX( I ) = ( DAT( 1,1)  +  DAT (1,2))*. 5 
BX ( N+ 1  ) =  BX ( N+1  )  +  BX ( I , J  ) 

DX (N  +  1  )  =  DX ( N  +  1  )  +  DX ( I , J  ) 

200  CONTINUE 
RETURN 
END 
C 
C 

C  SUBROUTINE  SOJOUP 
C 

SUBROUTINE  SOJOURf  SEED2 , A2 ,W,K) 
INTEGER  I , W , K 
REAL  A2(W,K) 

DOUBLE  PRECISION  EFF , SEED 2 
EFF  =  21 47483647. 0D0 
DO  10  1=1 ,K 

SEED2=DM0D( 16807. 0D0  *  SEED2 , EFF ) 
A2( 1  , I  )  =  SEED2/EFF 
A2 ( 2  ,  I )  «  1.0  -  A2( 1  ,1) 

10  CONTINUE 
RETURN 
END 
C 
C 

C  SUBROUTINE  STATE 
C 

SUBROUTINE  STATE( SEED1  , A1  , W , K ) 
INTEGER  I.K.W 
REAL  A1 (W.K) 

DOUBLE  PRECISION  EFF , SEED1 
EFF=  21 47483647. 0D0 
DO  10  1=1 ,K 

SEED1 =DMOD( 1 6807 . 0D0  *  SEED1 , EFF ) 
A 1  (  ,  I  )  =  SEED1  /EFF 

A 1 ( 2  ,  I )  =  1.0  -  A 1 ( 1  , I ) 

10  CONTINUE 
RETURN 
END 
C 
C 
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C  SUBROUTINE  STAT 
C 

SUBROUTINE  STAT ( N , R , BX , DX , BOX , DOX ) 

TMTTrmrp  TPM 

REAL  BX(N+1  )  ,DX(N+1  ) ,BOX(R) ,DOX(R) 

BOX ( 2 ) »  0.0 
DOX ( 2 ) =  0.0 
C 

C  COMPUTE  THE  ESTIMATES  OF  THE  SAMPLE  MEAN  AND 
C  VARIANCE 
C 

BOX ( 1 )  =  BX ( N+ 1 , J )/N 
DOX (  1  )  =  DX (  N  t- 1  ,  J  )/N 
DO  260  1  =  1  ,  N 

BOX ( 2  )  =  BOX ( 2  )  +  (BX(I)-BOX(I  ))*  *2 
DOX ( 2  )  =  DOX ( 2  )  +  (DX(I)-DOX(I  ))**2 
260  CONTINUE 

BOX ( 2  )  =  B0X(2)/(N*(N-1  )) 

DOX ( 2  )  =  DOX (2)/(N*(N-1  )) 

RETURN 

END 


M 
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APPENDIX  D 


FORTRAN  PROGRAM  LISTING  OF  THE  BCD  SIMULATION  USING 

STRATIFIED  SAMPLING 


DIMENSION  BX( 1 01  ) ,DX( 1 01  )  , SEED( 2 ) , BOX( 2  )  ,D0X(2) 

INTEGER  I,N,BB,DD,R 

REAL *4  X,Y,Z,BX,DX,TXT,BOX, DOX 

DOUBLE  PRECISION  SEED 

DATA  SEED  / 1 234 . 0D0 , 567890 1 23 . 0D0/ 

R*  2 
N-1  00 
C 

C  RECEIVE  INPUT  DATA  FROM  TERMINAL 
C 

WRITE( *  ,3) 

3  FORMAT ( IX, ’ENTER  THE  NUMBER  OF  BOMBERS’) 

READ  ’  (  12  )  ’  ,  BB 

WRITE( »,4) 

4  FORM.  .1 v  1 X ,  ’ ENTER  THE  RATE  WHICH  A  BOMBER  SHOOTS’ 
v*  '  DOWN  A  DEFENDER  ’  ) 

READ  ’ (F5.3) ’ ,  Y 
WRITE( * ,5 ) 

5  FORMAT ( IX, ’ENTER  THE  NUMBER  OF  DEFENDERS’) 

READ  ’  ( 12)  ’  ,  DD 

WRITE( * .6) 

6  FORMAT ( IX, ’ENTER  THE  RATE  WHICH  A  DEFENDER’ 

&  ’SHOOTS  DOWN  A  BOMBER’) 

READ  ’ ( F5 . 3  )  ’ ,  X 

WRITE( • ,7) ’ENTER  THE  RATE  WHICH  FREE  DEFENDERS’ 
&  ’FIND  FREE  BOMBERS’ 

7  FORMAT ( 1 X , A ) 

READ  ’ ( F5 . 3 ) ’ ,  Z 

WRITE( * ,8) ’ENTER  THE  TIME  DURATION  OF  THE’ 

&  ’BATTLE’ 

8  FORMAT ( 1 X , A ) 

READ  ’ (F5.3) ’ ,  TXT 
C 

C  RUN  REPLICATIONS  OF  N  BATTLES  AND  OBTAIN  SUMMARY  OF 
C  N  BATTLES 

C 

CALL  BATTLE( N , R , SEED , X , Y , Z , BB , DD , TXT , BX , DX ) 

C 

C  COMPUTE  STATISTICAL  OUTPUT  ANALYSIS  OBTAIN  PARAMETER 
C  ESTIMATES 
C 

CALL  STAT( N , R , BX , DX , BOX , DOX ) 


c 

C  FORMAT  AND  PRINT  OUTPUT  OF  PARAMETER  ESTIMATES 
C 

WRITE(3.279)  N 

279  FORMAT ( 1 5X, ’SAMPLE  SIZE’ , 16 , / 1 5X , 6( ), 

&  IX. 4(  ’-’  )) 

WRITE( 3 . 280 )  BB.DD 

280  &  FORMAT( //IX, 14, 2X, ’BOMBERS’ ,6X, ’VERSUS’ , 
&  6X , 14 , 2X , ’DEFENDERS ’ , 

&  / 1 X , 1 3 (  ) , 1 8X , 1 4 (  ’-’  )) 

290  WRITE  ( 3 , 300  )TXT 

300  FORMAT ( /  /  1 8X ,  ’TIME’  , F6 . 1  ,  /  ,  47(  ’ - ’  ) ) 

WR ITE( 3.310) 

3 1 0  FORMAT (  1 9X ,  ’ BOMBER ’  .  2X  ,  ’ DEFENDER ’  . 

&  / 1 8X , 7(  ’  * ’  ),2X,8(  ’  =  ’  )) 

WR ITE( 3,320)  BOX ( 1  ) , DOX ( 1  ) 

320  FORMAT ( IX. ’AVERAGE’ .7X.2F10.4) 

WRITE( 3 , 330 )  BOX ( 2 ) , DOX ( 2 ) 

330  FORMAT ( IX, ’VARIANCE’ .6X.2F10.4) 

STOP 

END 

C 

C 

C  SUBROUTINE  BATTLE 
C 

SUBROUTINE 

&  BATTLE( N , R , SEED , X , Y , Z , BB . DD , TXT , BX , DX ) 

INTEGER  BB , DD ,H,I,J,K,N,R,W 

REAL  GA( 3 , 1 000 ) , SO J (3,1 000 ) , BM( 34.3), 

&  DF( 34 , 3  )  , BX( N+1  ),DX(N+1  ), X , Y , Z , T , TXT , TIME , 

&  INF.NC.BW.DW.B.C.D.BM.DF 
TOUBLE  PRECISION  SEED(R) 

K  =  1000 

W  =  3 
C 

C  INITIALIZE  STATISTICAL  COUNTERS 
C 

BX ( N+ 1  )=  0.0 
DX (N+1  )»  0.0 

C 

C  RUN  N  REPLICATIONS 
C 

DO  200  1-1 ,N 

CALL  SOJOUR( SEED( 1 ) , SOJ , W, K ) 

CALL  STATE ( SEED( 2 ) , GA , W , K ) 

DO  175  J- 1 ,W 
C 

C  INITIALIZE  START-TO-BATTLE  VALUES 
C 


H  =  0 
T  =  0 . 0 


97 


OO  OOO  GOO  GOG  GGGG  GGG 


C  -0.0 

B-REAL(BB) 

D-REAL(DD) 

NC-  1 .  0 
BW-  0.0 
DW=  0 . 0 

OBSERVE  AN  OCCURRENCE  OF  AN  EVENT 
100  H-H+1 

DETERMINE  NEXT  INTERIM  EVENT  AND  UPDATE  THE  STATE 
(B,C,D) 

IF( GA( J , H )  .LE.  BW)  THEN 
B-B  +  1.0 
C-C  -  1.0 
D-D 

ELSE  IF  ( GA (  J  ,  H  )  . LE .  ( BW+NC ) )  THEN 
B-B  -1.0 
C=C  +1.0 
D=D  -  1.0 
ELSE 
B»B 

C-C  -  1.0 
D-D+  1 . 0 
END  IF 

COMPUTE  MEAN  TIME  IN  STATE  (B,C,D) 

IF( ( B  .EQ.  0.0  .OR.  D  .EQ.  0.0). AND.  C 
&  .EQ.  0.0)  THEN 

INF*  1000000.0 
ELSE 

INF-  1  . 0/ ( Z  *B*D  +  (X+Y)*C) 

END  IF 

COMPUTE  SOJOURN  TIME  OF  STATE  (B,C,D) 

TIME-  -INF  *  ALOG( SO J ( J , H ) ) 

ADVANCE  SIMULATED  TIME  OF  THE  AIR  BATTLE 
T-  T  +  TIME 

COMPUTE  PROBABILITY  OF  NEXT  INTERIM  EVENTS  OCCURRING 
C 

NC=  Z  *  B*D* INF 
BW-  Y*C  *  INF 
DW-  X*C»INF 
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c 

C  CHECK  FOR  OCCURRENCE  OF  TERMINATION  EVENT 
C 

IF  (T  .LT.  TXT(TI))  GOTO  100 
C 

C  RECORD  RESULTS  OF  BATTLE 
C 

BM(  I ,  J  )  =  B  +  C 
DF( I , J  )«  D  +  C 
175  CONTINUE 

C 

C  ACCUMULATE  SUMMARY  OF  N  BATTLE  OUTCOMES 
C 

BX (  I  )  =  ( BM( 1,1)  +  BM( 1,2)  +  BM( I , 3  )  / 3 . 0 
DX (  I  )  =  (DF(I.I)  +  DF (1,2)  +  BM(I,3)  . .0 
BX (N  +  1  )  =  BX ( N  +  1  )  +  BX ( I ) 

DX ( N+ 1  )*  DX ( N+ 1  )  +  DX ( I ) 

200  CONTINUE 
RETURN 
END 
C 
C 

C  SUBROUTINE  SOJOUR 
C 

SUBROUTINE  SO JOUR( SEED2 , B2 , W , K ) 

INTEGER  I , W , K , J 
REAL  B2(W,K),A2 
DOUBLE  PRECISION  EFF.SEED2 
EFF*  214748364-7. 0D0 
DO  10  1=1 , K 

SEED2=DM0D( 1 6807.0D0  *  SEED2 , EFF ) 

A2=  SEED2/EFF 
DO  5  J= 1 ,W 

B2( J , I ) =  AMOD(A2  +  ((J-1)  *  1  . 0  ) / 3  -  0 , 1  . 0 ) 
5  CONTINUE 

10  CONTINUE 
RETURN 
END 
C 
C 

C  SUBROUTINE  STATE 
C 

SUBROUTINE  STATE( SEED1 , A1 , W , K ) 

INTEGER  I , K , W , J 
REAL  A1(W,K),A2 
DOUBLE  PRECISION  EFF , SEED1 
EFF=  2147483647. 0D0 
DO  10  I-1 ,K 

SEED1=DMOD( 16807. 0D0  *  SEED1 , EFF ) 

A2=  SEED1/EFF 
DO  5  J= 1 ,W 
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A1 ( J , I )  ■  AM0D(A2  +  (CJ-1)  »  1 .0)/3.0, 1 .0) 
CONTINUE 
RETURN 
END 


SUBROUTINE  STAT 

SUBROUTINE  STAT ( N , R , BX , DX , BOX , DOX ) 

INTEGER  R  N 

REAL  BX( N+1 ) , DX( N+1 ) ,BOX(R) ,DOX(R) 

BOX ( 2  )  ■  0.0 
DOX  (  2 ) •  0.0 

COMPUTE  THE  ESTIMATES  OF  THE  SAMPLE  MEAN  AND 
VARIANCE 

BOX ( 1 )  «  BX( N+1 )/N 
DOX ( 1 )  =  DX( N+1 )/N 
DO  260  I- 1  ,N 

BOX ( 2)-  BOX ( 2 )  +  (BX(I)-BOX(I ))*»2 
DOX(2)»  DOX ( 2  )  +  (DX(I)-DOX(I  ))**2 
260  CONTINUE 

BOX ( 2  )  =  BOX ( 2)/(N»(N-1  )) 

DOX ( 2  )  ■  DOX ( 2)/(N*(N-1  )) 

RETURN 

END 
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STATISTICAL  OUTPUT  DATA  FROM  THE  BCD  SIMULATIONS 


TABLE  E.1  OUTPUT  PARAMETERS  OF  THE  BCD  MODEL  ESTIMATED 
FROM  CRUDE  SIMULATION 


SIMULATION 

SCENARIO 

DEFENDER 

BOMBER 

DEFENDER 

BOMBER 

TIME 

MEAN 

VAR 

MEAN 

VAR 

1 

1  0 

1  0 

25 

7 . 6 

.0191 

7.5 

.  0250 

2 

1  0 

1  0 

75 

4 . 4 

.0387 

4.5 

.0358 

3 

1  0 

1  0 

1  25 

2.8 

.0375 

3.4 

.0417 

4 

1  0 

30 

25 

5.3 

.0328 

25 . 1 

.  0541 

5 

1  0 

30 

75 

1 .3 

.0157 

21.3 

.  1  402 

6 

1  0 

30 

1  25 

.  3 

.0036 

20.8 

.  1  725 

7 

1  0 

50 

25 

4 . 2 

.0251 

44.7 

.  0558 

8 

1  0 

50 

75 

.6 

.  0064 

40.8 

.  1  654 

9 

1  0 

50 

1  25 

.  1 

.0014 

39.7 

.1874 

1  0 

30 

1  0 

25 

25. 1 

.  0458 

5.4 

.0350 

1  1 

30 

1  0 

75 

21.1 

.  1  407 

1  .  4 

.0141 

1  2 

30 

1  0 

125 

19.8 

.2247 

.  4 

.  0045 

1  3 

30 

30 

25 

18.7 

.0933 

18.5 

.0842 

1  4 

30 

30 

75 

8.2 

.1618 

8.4 

.  1895 

15 

30 

30 

125 

".7 

.  1  984 

4.7 

.  1  341 

1  6 

30 

50 

25 

14.4 

.0628 

34.5 

.1104 

17 

30 

50 

75 

3.4 

.0362 

22 . 9 

.3479 

18 

30 

50 

1  25 

.  9 

.0136 

20 . 1 

.  4881 

19 

50 

1  0 

25 

44.4 

.0651 

4 . 1 

.  0317 

20 

50 

1  0 

75 

40.7 

.  1  427 

.  7 

.  0063 

21 

50 

1  0 

125 

41  .  1 

.  1  666 

.  1 

.0014 

22 

50 

30 

25 

34 . 3 

.  1  687 

14.0 

.  0826 

23 

50 

30 

75 

23.5 

.3339 

3.1 

.0359 

24 

50 

30 

125 

21  .  1 

.5300 

.8 

.0107 

25 

50 

50 

25 

27 . 9 

.  1  590 

26.4 

.1501 

26 

50 

50 

75 

9.5 

.2717 

11.5 

.2865 

27 

50 

50 

1  25 

5.2 

.2126 

7.9 

.  2986 
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TABLE  E.2  OUTPUT  PARAMETERS  OF  THE  BCD  MODEL  ESTIMATED 
FROM  THE  SIMULATION  USING  ANTITHETIC  VARIATES 


SIMULATION 

SCENARIO 

DEFENDER  BOMBER 

TIME 

DEFENDER 

MEAN  VAR 

BOMBER 

MEAN  VAR 

1 

10 

10 

25 

7.6 

.  0067 

7 . 5 

.0105 

2 

10 

10 

75 

4 . 5 

.  0048 

4 . 4 

.  0060 

3 

10 

10 

125 

3.2 

.  0079 

3.3 

.  0057 

4 

10 

30 

25 

5 . 1 

.  0129 

25 . 1 

.  0100 

5 

10 

30 

75 

1 . 2 

.  0071 

21 . 3 

.  0329 

6 

10 

30 

125 

.  3 

.0024 

20.2 

.  0446 

7 

10 

50 

25 

4.4 

.  0123 

44 . 4 

.0146 

8 

10 

50 

75 

.  6 

.  0060 

40.9 

.  0395 

9 

10 

50 

125 

.  1 

.0013 

40 . 1 

.  0683 

10 

30 

10 

25 

25.4 

.  0116 

5 . 0 

.0116 

11 

30 

10 

75 

21 . 1 

.0271 

1.2 

.  0056 

12 

30 

10 

125 

20.2 

.  0558 

.  3 

.  0030 

13 

30 

30 

25 

18.5 

.0164 

18.4 

.  0210 

14 

30 

30 

75 

8.5 

.0251 

8.4 

.0270 

15 

30 

30 

125 

5.1 

.0307 

5 . 1 

.0352 

16 

30 

50 

25 

14 . 6 

.  0223 

34.6 

.0308 

17 

30 

50 

75 

3.0 

.0138 

23.6 

.  0542 

18 

30 

50 

125 

1 . 0 

.0075 

20.6 

.0446 

19 

50 

10 

25 

44.0 

.0169 

4 . 3 

.0104 

20 

50 

10 

75 

40.2 

.0410 

.6 

.0059 

21 

50 

10 

125 

40.6 

.0451 

.  1 

.0015 

22 

50 

30 

25 

34.3 

.0279 

14.4 

.0368 

23 

50 

30 

75 

23.1 

.0375 

3.1 

.0194 

24 

50 

30 

125 

20.5 

.0739 

.  9 

.  0077 

25 

50 

50 

25 

27.3 

.  0266 

27.1 

.0357 

26 

50 

50 

75 

10.5 

.0269 

10.7 

.0309 

27 

50 

50 

125 

6.4 

.0458 

6 . 4 

.  0472 
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TABLE  E . 3  EFFICIENCY  OF  ANTITHETIC  VARIATES  IN  THE  BCD 
MODEL 


SIMULATION 

SCENARIO 

DEFENDER 

BOMBER 

DEFENDER 

BOMBER 

TIME 

RE 

IP(*) 

RE 

IP(*) 

1 

10 

10 

25 

2 . 9 

65 . 0 

2.4 

58.0 

2 

10 

10 

75 

8.1 

87 . 6 

6 . 0 

83.2 

3 

10 

10 

125 

4 .7 

78.9 

7.3 

86.3 

4 

10 

30 

25 

2 . 5 

60 . 7 

5 . 4 

81 . 5 

5 

10 

30 

75 

2 . 2 

54.8 

4 . 3 

76.5 

6 

10 

30 

125 

1 . 5 

33.3 

3.9 

74 . 1 

7 

10 

50 

25 

2 . 0 

51 . 0 

3.8 

73.8 

8 

10 

50 

75 

1  .  1 

6.3 

4.2 

76. 1 

9 

10 

50 

125 

1  .  1 

7. 1 

2.7 

63 . 6 

10 

30 

10 

25 

3.9 

74 . 7 

3.0 

66.9 

11 

30 

10 

75 

5 . 2 

80.7 

2.5 

60.3 

12 

30 

10 

125 

4 . 0 

75 . 2 

1 . 5 

33.3 

13 

30 

30 

25 

5.7 

82.4 

4.0 

75 . 1 

14 

30 

30 

75 

6.4 

84 . 5 

7.0 

85.8 

15 

30 

30 

125 

6.5 

84 . 5 

3.8 

73.8 

16 

30 

50 

25 

2.8 

64 . 5 

3.6 

72. 1 

17 

30 

50 

75 

2.6 

61 . 9 

6.4 

84 . 4 

18 

30 

50 

125 

1.8 

44 . 9 

10.9 

90.9 

19 

50 

10 

25 

3.9 

74 . 0 

3.0 

67.2 

20 

50 

10 

75 

3.5 

71.2 

1 . 1 

6.3 

21 

50 

10 

125 

3.7 

72 . 9 

.9 

7.1 

22 

50 

30 

25 

6.0 

83.4 

2.2 

55.4 

23 

50 

30 

75 

8.9 

88.8 

1.9 

46.0 

24 

50 

30 

125 

7 . 2 

86.1 

1.4 

28.0 

25 

50 

50 

25 

6.0 

83.3 

4.2 

76.2 

26 

50 

50 

75 

10.1 

90.1 

9 . 3 

89.2 

27 

50 

50 

125 

4.6 

78.5 

6 . 3 

84.2 
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TABLE  E.4 

OUTPUT  PARAMETERS 

OF  1 

THE  BCD 

MODEL 

ESTIMATED 

FROM  SIMULATION  USING 

STRATIFIED 

SAMPLING 

SCENARIO 

DEFENDER 

BOMBER 

SIMULATION 

DEFENDER 

BOMBER 

TIME 

MEAN 

VAR 

MEAN 

VAR 

iB 

1 

1  0 

1  0 

25 

7.4 

.0119 

7.5 

.0114 

Jt1*! 

2 

1  0 

1  0 

75 

4.4 

.0248 

4.5 

.  0237 

3 

1  0 

1  0 

125 

3.1 

.0160 

3.2 

.0197 

1 

4 

1  0 

30 

25 

5.1 

.0124 

25.0 

.0252 

5 

1  0 

30 

75 

1  .2 

.0076 

21  .2 

.0710 

6 

1  0 

30 

125 

.4 

.0038 

20.5 

.  0801 

% 

7 

1  0 

50 

25 

4.4 

.0126 

44 . 0 

.0375 

8 

1  0 

50 

75 

.9 

.0083 

40 . 2 

.1021 

9 

1  0 

50 

125 

.  1 

.  0006 

40 . 0 

.  0930 

10 

30 

1  0 

25 

25 . 1 

.0233 

5.3 

.0124 

1  1 

30 

10 

75 

21.1 

.  0726 

1  .3 

.0118 

1  12 

30 

1  0 

125 

20.3 

.  0877 

.4 

.  0027 

s 

13 

30 

30 

25 

18.4 

.0283 

18.4 

.0193 

14 

30 

30 

75 

8.0 

.0439 

8.2 

.  0700 

!/ 

15 

30 

30 

125 

5.5 

.0505 

4.8 

.  0426 

L 

16 

30 

50 

25 

14.2 

.0219 

34.5 

.0798 

'  17 

30 

50 

75 

3.2 

.021  1 

23.3 

.  1  236 

v: 
’  ► 

18 

30 

50 

125 

.8 

.  0081 

20. 1 

.  1880 

19 

50 

10 

25 

44.6 

.0339 

4.3 

.0163 

& 

20 

50 

1  0 

75 

41  .0 

.0872 

.7 

.0053 

& 

21 

50 

10 

125 

40.4 

.  1  089 

.  1 

.0012 

1  22 

50 

30 

25 

34.6 

.0571 

14.3 

.0439 

23 

50 

30 

75 

22.6 

.1187 

3.3 

.0197 

;<J 

i  24 

50 

30 

125 

20.9 

.2402 

1  .0 

.0086 

v; 

1  25 

50 

50 

25 

27.7 

.  0630 

27.0 

.  1  084 

i  26 

50 

50 

75 

10.5 

.  0898 

10.7 

.  0881 

•S 

i  27 

i 

50 

50 

125 

6.5 

.  091  1 

6.4 

.  1  084 

TABLE  E . 5  EFFICIENCY  OF  STRATIFIED  SAMPLING  IN  THE  BCD 
MODEL 


SIMULATION 

SCENARIO 

DEFENDER 

BOMBER 

DEFENDER 

BOMBER 

TIME 

RE 

IP(*) 

RE 

IP(%) 

1 

1  0 

1  0 

25 

1  .6 

39.0 

2.2 

54.4 

2 

1  0 

1  0 

75 

1  .  6 

36 . 0 

1  .5 

33.8 

3 

10 

1  0 

125 

2.3 

57.3 

2 . 1 

52.8 

4 

1  0 

30 

25 

2.6 

62.2 

2 . 1 

53.4 

5 

1  0 

30 

75 

2 . 1 

51  .  6 

2 . 0 

49.4 

6 

1  0 

30 

125 

.9 

-5.5 

2.2 

53 . 6 

7 

1  0 

50 

25 

2.0 

49.8 

1  .5 

32.8 

8 

1  0 

50 

75 

.8 

-29.7 

1  .6 

38.3 

9 

1  0 

50 

1  25 

2.3 

57. 1 

2.0 

50 . 4 

1  0 

30 

1  0 

25 

2.0 

49 . 1 

2.8 

64.6 

1  1 

30 

1  0 

75 

1  .9 

48.4 

1  .  2 

16.3 

1  2 

30 

1  0 

125 

2.6 

61  . 0 

1  .7 

40 . 0 

1  3 

30 

25 

3.3 

69.7 

4.4 

77. 1 

1  4 

30 

30 

75 

3.7 

72.9 

2.7 

63 . 1 

1  5 

30 

✓  — 

125 

3.9 

74 . 5 

3 . 1 

68.2 

1  6 

30 

50 

25 

2.9 

65 . 1 

1  .4 

27.7 

1  7 

30 

50 

75 

1  .7 

41.7 

2.8 

64.5 

18 

30 

50 

125 

1  .7 

40.7 

2 . 6 

61  .5 

19 

50 

10 

25 

1  .9 

47.9 

1  .9 

48.6 

20 

50 

1  0 

75 

1  .6 

39.0 

1  .2 

15.9 

21 

50 

10 

125 

1  .5 

34.6 

1  .2 

14.2 

22 

50 

30 

25 

3.0 

66.2 

1  .9 

46.9 

23 

50 

30 

75 

2.8 

64.5 

1  .8 

45 . 1 

24 

50 

30 

125 

2.2 

54.7 

1  .2 

19.6 

25 

50 

50 

25 

2.5 

60.4 

1  .4 

27.8 

26 

50 

50 

75 

3.0 

66.9 

3.3 

69 . 2 

27 

50 

50 

125 

2.3 

57.1 

2.8 

63 . 7 

1  05 


LIST  OF  REFERENCES 


Al-Zayani,  Abdul-Latif  Rashid,  Formulation  and  Analysis 
Combat-Logistics  Problems,  Ph.d.  Thesis,  Naval 
Postgraduate  School,  Monterey,  California,  1986. 

Andrtasson,  Ingar  J.,  "Antithetic  Methods  in  Queuing 
Simulations,"  Report  NA.  72.58,  Department  of 
Computer  Science,  Royal  Institute  of  Technology, 
Stockholm, 1 972. 

Bratley,  Paul;  Fox,  Bennett  L.:  and  Schrage,  Linus  E.; 

A  Guide  to  Simulation ,  Second  Edition,  Springer- 
Verlag,  New  York,  1987. 

Cheng,  R.  C.  H.,  "Variance  Reduction  Methods," 
Proceedings  of  the  1986  Vinter  Conference. 

Institute  of  Electrical  and  Electronic  Engineers, 
Washington,  D.  C.,  pp .  60-68. 

Fishman,  George  S. ,  Principles  of  Discrete  Event 
Simulation ,  Wiley,  New  York,  1978. 

Gaver ,  Donald  P.  and  Thompson,  Gerald  L.,  Programming 
and  Probability  Models  In  Operations  Research, 
Brooks/Cole  Publ ishing  Company ,  Monterey ,  ~ 

California,  1973. 

Hammersley,  J.  M.  and  Handscomb ,  D.  C.,  Monte  Carlo 
Methods .  Methuen  &  Co  Ltd,  London,  1964. 

Handscomb,  D.  C.,  "Monte  Carlo  Techniques: 

Theoretical,"  In  The  Design  of  Computer  Simulation 
Experiments ,  pp .  252-262.  Edited  by  Thomas  H. 

Naylor,  Duke  University  Press,  Durham,  N.C.,  1969. 

Hartman,  James  K.  "Lectures  Notes  in  High  Resolution 
Combat  Modelling,"  Naval  Postgraduate  School, 
Monterey,  California,  1985. 

John,  Peter  W.  M. ,  Statistical  Design  and  Analysis  of 
Experiments ,  The  Macmillan  Company,  New  York,  1971. 

Johnson,  Thomas  C.;  Bates,  Carl  B.;  and  Graham,  Breton 
C . ,  Application  of  Antithetic  Variates  to  the  CQSAGE 
IV  Mode  1 ,  CAA-TP-85-4,  U.  S.  Army  Concepts  Analysis 
Agency,  Bethesda,  Maryland,  1985. 


1  06 


MiMMMiiiiiiMiaMiMiaiiH 


Kleijnen,  Jack  P.  C.,  Statistical  Techniques  In 

Slmulat ion .  Part  I ,  Marcel  Dekker,  New  York,  1974. 

Larson,  Harold  J.,  Introduction  to  Probability  Theory 
and  Statistical  Inference.  John  Wiley  &  Sons,  New 
York,  1982. 


Law,  Averill  M.  and  Kelton,  David,  Simulation  Modeling 
and  Analysis.  McGraw-Hill,  New  York,  1982. 

Lewis,  P.  A.  W.  and  Orav  E.  J.  "Simulation  Methodology 
for  Statisticians,  Engineers,  and  Operations 
Analysts."  Naval  Postgraduate  School,  Monterey, 
California,  1985. 

McGrath,  E.  J.  and  Irving,  D.  C.,  "Application  of 
Variance  Reductionto  Large  Scale  Simulation 
Problems,"  Compute,  and  Ops.  Res.,  Vol  I,  Pergamon 
Press,  Oxford,  1974,  pp .  283-311. 

Morgan,  Byron  J.  T.  ,  Elements  of  Simulation.  Chapman 
and  Hall,  New  York,  1986. 

Moy .  William  A.,  "Monte  Carlo  Techniques:  Practical," 
In  The  Design  of  Computer  Simulation  Experiments, 
pp .  263-288.  Edited  by  Thomas  H.  Naylor,  Duke 
University  Press,  Durham,  N.  C.,  1969. 

Nelson,  Barry  L.,  "A  Decomposition  Approach  to  Variance 
Reduction,"  Proceedings  of  the  1985  Winter 
Conference ,  Institute  of  Electrical  and  Electronic 
Engineers,  Washington,  D.  C.,  pp .  23-31. 


1  07 


INITIAL  DISTRIBUTION  LIST 


No. 

1.  Defense  Technical  Information  Center 
Cameron  Station 

Alexandria,  Virginia  22304-6145 

2.  Library,  Code  0142 
Naval  Postgraduate  School 
Monterey,  California  93943-5002 

3.  Director 

U.  S.  Army  Concepts  Analysis  Agency 
ATTN:  Mr.  Carl  Bates 
8120  Woodmont  Avenue 
Bethesda,  Maryland  20814-2797 

4.  Professor  Donald  P.  Gaver,  Code  55Gv 
Department  of  Operations  Research 
Naval  Postgraduate  School 
Monterey,  California  93943-5000 

5.  Professor  Alvin  F.  Andrus,  Code  55As 
.Department  of  Operations  Research 

Naval  Postgraduate  School 
Monterey,  California  93943-5000 

6.  Professor  Sam  H.  Parry,  Code  55Py 
Department  of  Operations  Research 
Naval  Postgraduate  School 
Monterey,  California  93943-5000 

7.  Director 

LSA  TRAC-FLVN 

ATRC-FSC  ( CPT  Curtis  Smith) 

FT  Leavenworth,  Kansas  66027-5220 


Copie 

2 


2 


3 


2 


1 


1 


4 


108 


